{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import the basic libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import torch\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython import display\n",
    "import os\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize the environment for running the experiment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.makedirs('./results/HGaussian_MINEE/', exist_ok=True)\n",
    "\n",
    "# Paths for loading/storing results.\n",
    "name = 'results/HGaussian_MINEE' # filename\n",
    "chkpt_name = name+'.pt'              # checkpoint\n",
    "fig_name = name+'.pdf'               # output figure\n",
    "\n",
    "# use GPU if available\n",
    "if torch.cuda.is_available(): \n",
    "    torch.set_default_tensor_type(torch.cuda.FloatTensor)\n",
    "else:\n",
    "    torch.set_default_tensor_type(torch.FloatTensor)\n",
    "\n",
    "# initialize random seed\n",
    "np.random.seed(0)\n",
    "torch.manual_seed(0);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate data using the gaussian model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_size = 400   # sample size\n",
    "rho = 0.9           # model parameter\n",
    "p = 0.8             # parameter of bernoulli distribution\n",
    "mu = 1              # mean of Gaussian distribution\n",
    "rep = 1             # number of repeated runs\n",
    "d = 1               # number of dimensions for X (and Y)\n",
    "X = np.zeros((rep,sample_size,d))\n",
    "Y = np.zeros((rep,sample_size,d))\n",
    "Z = np.zeros((rep,sample_size,d))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(rep):\n",
    "    Z[i] = np.random.normal(\n",
    "            loc=mu,\n",
    "            scale=rho,\n",
    "            size=sample_size).reshape(sample_size,1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate data from bernoulli distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import bernoulli\n",
    "for i in range(rep):\n",
    "    X_ = bernoulli.rvs(p, size=sample_size)\n",
    "    X_[np.where(X_==0)] = -1\n",
    "    X[i] = X_.reshape(sample_size, d)\n",
    "\n",
    "# data_Y is the product of data from bernoulli distribution and mixed Gaussian distribution.\n",
    "for i in range(rep):\n",
    "    for j in range(sample_size):\n",
    "        Y[i,j] = X[i,j] * Z[i,j]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A plot of the first dimension of $Y$ against that of $X$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEWCAYAAABv+EDhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAYbklEQVR4nO3de7RkVX3g8e+vebQPQGjpCCiXBkQMw7iEuYM4SUSRsUGNGESDMYrx0QtNXGPGgI9GYwQTHrOIa0ZH6UwI4gMfCJHYkhZE4nLSMGlciA2INCgPeXQrtsgjDU3/5o86x65bt+6zq+pU9f5+1qrVtc9r/+6p6vqdvc85+0RmIkkqz4KmA5AkNcMEIEmFMgFIUqFMAJJUKBOAJBXKBCBJhTIBaOAi4pqIeMeA6npXRDwQEQ9HxDNnsfxPI+KYQcQ2DCLiwog4s+k41AwTgPqi+iF9rPrhfSAi/iEidpnjNpZEREbEjvOMYSfgPOAVmblLZv5iPtuZZvsZEc/t5TalQTIBqJ9+PzN3AQ4H/jNw+oDrfxbwFOCmAdcrjQQTgPouM38GXAEc2jkvIhZExOkRcWdErI+IiyLiGdXs71b/bqxaEi/usv7CiPhERNxbvT5RTXsecGvb+ld3iy0i3lzV/YuIWN4x74iIWB0RGyPivoj4ZETsXM2rY/tBFdsfRsQeEfGNiNgQEb+s3j9nqv0SEe+PiJ9FxK8j4taIePlM9VbzMyLeHRG3VeueEREHVus8FBFfaYvzpRFxT0R8KCJ+XrXM3jRNTK+OiBuquv81Il4wU7waYZnpy1fPX8BPgWOq9/vSOgo/oypfA7yjev82YB1wALALcCnwuWreEiCBHaep52PAtcBvAYuBf22rZ9r1gUOAh4GXAAtpdRdtbov7PwFHAjtW27oFeG/b+gk8t638TOB1wNOAXYGvAv84Rd0HA3cD+7TFeuAc6r0c2A34D8Am4NvVPnwGcDNwcrXsS6u/6bzqbzwKeAQ4uJp/IXBm9f5wYD3wImAH4OTqc1w4Xby+RvdlC0D99I8RsRH4HvAvwF93WeZNwHmZeUdmPgx8EDhpDv3+bwI+lpnrM3MD8FfAm2e57onANzLzu5m5CfgwsKWemZnXZ+a1mbk5M38KnE/rB7SrzPxFZn4tMx/NzF8DH59m+Sdp/bAeEhE7ZeZPM/P2OdR7dmY+lJk3AWuBb1X78Fe0WluHdSz/4czclJn/AqwE3tAlpncC52fmdZn5ZGZ+llZyOXK6eDW6TADqp9dm5u6ZuV9mvjszH+uyzD7AnW3lO2kd+T5rlnV0W3+fOax7d13IzEeA35wojojnVd0490fEQ7QS2J5TbSwinhYR51ddSg/R6sLaPSJ26Fw2M9cB7wU+CqyPiC9FxD5zqPeBtvePdSm3n3D/ZfW31abaR/sB76u6fzZWyXtfWkf9U8ar0WUCUNPupfXDUxuj1WXxAK2ujvmsf+8s676P1g8c0PoBp9WNU/s08CPgoMzcDfgQENNs7320ukpeVC3/knrT3RbOzC9m5u9W8Sdw9jzrnckeEfH0tvJU++hu4ONV0q5fT8vMi2eIVyPKBKCmXQz8eUTsX10m+tfAlzNzM7CBVpfMATOsf3pELI6IPYGPAJ+fZd2XAK+OiN+tTpp+jIn/J3YFHgIejojnA+/qWP+Bjth2pXX0vTEiFgF/OVXFEXFwRBwdEQuBf6/We3KW9c7HX0XEzhHxe8CraZ2f6PR3wCkR8aJoeXpEvCoidp0hXo0oE4CadgHwOVrdJT+h9ePyHoDMfJRWP/r/rbokjuyy/pnAGuBG4IfA96tpM6r6z/8U+CKt1sAvgXvaFvkL4I+AX9P6cfxyxyY+Cny2iu0NwCeApwI/p3Vi+p+nqX4hcFa17P20TmJ/aJb1ztX9tP62e4EvAKdk5o86F8rMNbTOA3yyWn4d8NZZxKsRFZk+EEbaXkXES4HPZ+aUl6OqXLYAJKlQJgBJKpRdQJJUKFsAklSoeY2y2JQ999wzlyxZ0nQYkjRSrr/++p9n5uLO6SOVAJYsWcKaNWuaDkOSRkpE3Nltul1AklQoE4AkFcoEIEmFMgFIUqFMAJJUKBOAJBXKBCBJQ+rUi1Zz6kWr+7Z9E4AkFWqkbgSTpBLUR/033vnghPK5b3lxT+uxBSBJhWqsBRART6H1FKiFVRyXZOaUj9CTpFLUR/rHnblyQrnXmuwC2gQcnZkPR8ROwPci4orMvLbBmCSpcSecswqALTmxfOlpS3taT2MJIFsPIni4Ku5UvXw4gSQNSKMngSNiB+B64LnApzLzui7LLAOWAYyNjQ02QElqQH2kX3cB9frIvzYUTwSLiN2By4D3ZObaqZYbHx9Ph4OWtL3rvAroBfstAuZ/LiAirs/M8c7pQ3EVUGZuBK4Bjm04FEkqRpNXAS0GnsjMjRHxVOAY4Oym4pGkYVEf6ffr+v9ak+cA9gY+W50HWAB8JTO/0WA8klSUJq8CuhE4rKn6JWnY3X7/Q33dvkNBSNKQemTT5r5u3wQgSUOmvvGrs7zd3AgmSequ88i/Xy2BobgMVJI0eCYASSqUXUCSNGTqO3877wTuNROAJA2Zzss/+3U5qF1AklQoWwCSNGQO3Gs3YGsXUF3uNROAJA2ZtXc9OG25V0wAkjRknrpz66e5vv6/LveaCUCShkx9x+/SM/r7QBhPAktSoWwBSNKQcSwgSSrUoMYCMgFI0pBZEK1/t+TEcq+ZACRpyBw6NnEoiLrcayYASRoy9TOA66uAtsdnAkuSuqgfBt9Z7nUiMAFI0pDxTmBJKpR3AktSoQZ1J7AJQJKGTP3D31le9eFX9bSexoaCiIh9I+I7EXFLRNwUEf+tqVgkqURNtgA2A+/LzO9HxK7A9RFxZWbe3GBMklSMxloAmXlfZn6/ev9r4Bbg2U3FI0mlGYrRQCNiCXAYcF2XecsiYk1ErNmwYcOgQ5Ok7VbjCSAidgG+Brw3Myc9+TgzV2TmeGaOL168ePABStKALYiJ4/90lnul0auAImInWj/+X8jMS5uMRZKGxXZ/H0BEBPD3wC2ZeV5TcUjSsBnUcNBNdgH9DvBm4OiIuKF6vbLBeCSpKI21ADLze0CfRrmWJM2k8ZPAkqRmmAAkqVAmAEkqlAlAkgplApCkQpkAJKlQJgBJKpQJQJIKZQKQpEKZACSpUEUlgFMvWs2pF61uOgxJGgpFJQBJ0laNPg9gUOqj/hvvfHBC+dy3vLixmCSpabYAJKlQRbQA6iN9j/wlaStbAJJUqCJaADWP/CVpK1sAklQoE4AkFaqoLiBJGgULqqelb8mJ5V4zAUjSkKl/+Kcq94pdQJJUqEYTQERcEBHrI2Jtk3FIUomabgFcCBw7qMocDE6Stmo0AWTmd4EHm4xBkkpVxElgB4OTpMma7gKaUUQsi4g1EbFmw4YNTYcjSduNoW8BZOYKYAXA+Pj4vC6GcjA4SZps6BNAL91+/0NNhyBJQ6Ppy0AvBlYDB0fEPRHx9n7Wd+Beu3HgXrv1swpJGhmNtgAy842DqMeTwJI02dCfBJYk9UcR5wA8CSxJk9kCkKRCFdECqHnkL0lb2QKQpEKZACSpUEUlgBPOWcUJ56xqOgxJGgpFJQBJ0lZFnASuj/of2bR5QvnS05Y2FpMkNc0WgCQVqogWQH2kv/SMlRPKkjSMnr6w9dNc91rU5V4rIgF0PgbSO4IlDbP6h3+qcq8UkQA6h4F2WGhJ8hyAJBWriBZA/QyAejhonwkgSbYAJKlYRbQAPAcgSZMVkQDsApKkyYpIAPXlnvUdwF7+KUmeA5CkYpkAJKlQRXQB1Xf+1nfTeSewJBWSALwKSJImm7ILKCK+GRFL+ll5RBwbEbdGxLqI+EA/65IkTTTdOYALgW9FxPKI2KnXFUfEDsCngOOAQ4A3RsQhva4H4LHHN/PY45unLEtSiabsAsrMr0TESuAjwJqI+BywpW3+edtY9xHAusy8AyAivgQcD9y8jduVJM3CTOcAngAeARYCu9KWAHrg2cDdbeV7gBd1LhQRy4BlAGNjY/Oq6NCxRcDWG8HqsiSVbMoEEBHHAucBlwOHZ+ajPa47ukzLSRMyVwArAMbHxyfNlyTNz3QtgOXA6zPzpj7VfQ+wb1v5OcC9fapLktRhunMAv9fnuv8NOCgi9gd+BpwE/FE/Kqqv9/f6f0mjYLt/JGRmbo6IPwNWATsAF/SxtSFJ6tDojWCZ+U3gm4OqzyN/SaNgUM8ELmosoFMvWj3pAfGSVKqiEoAkaasixgKqj/rr+wA8GSxpmC2oLpLfkhPLvVZEApCkUbIlpy/3ShEJoD7yn6osSSUqIgFI0ijZ7u8DkCR11zlacb9GL/YqIEkqVBEtgEE1pySpFzwJ3EMH7rUbsPXkb12WpJIVkQB8JrCkUeJ9AD1kC0DSKBnUQ6yKSAC2ACRpMq8CkqRCFdECuPS0pQAsPWPlhLIklayIBNA5BLSDwUlSIQnAcwCSRkl9cHrcmSsnlHutiAQwqNuqJakX6l6K+jLQfvVaFJEAJGmUDKrXoogE8NSdJw4FUZclaRgN6t4lfwklach4DkCSCjWocwDeCCZJhWqkBRARrwc+Cvw2cERmrulnfV4FJGmU1Ef6/b5nqakuoLXACcD5DdUvSUOv3/csNZIAMvMWgIg+jXHawauAJI2ifo9cPPS/hBGxDFgGMDY2Nq9t2AUkaZTUXT/1ZaAjdyNYRFwF7NVl1vLM/Ppst5OZK4AVAOPj4/N6MNqgHq8mSaOkbwkgM4/p17bnymcCSxol2/tJ4IHyiWCSNFlTl4H+AfC/gMXAyoi4ITP7Nki/o4FKGkX9HrK+qauALgMuG1R9tgAkabIiuoAGNa6GJPVSv88BOBSEJBWqiBbAoAZWkqReGPn7AIaJJ4ElabIiEkB9/f9UZUkaJt4H0EMLqiGH6i6gBYMZgkiShlpkjs64COPj47lmzfxHjl56RusqoFUfflWvQpKkoRcR12fmeOd0rwKSpCF16kWrf9MN1A9FdAGdcM6qruVLT+vbzceSNPSKSAAOBy1plHgZaA8dOrYI2Loz67IklayIBCBJo8TLQHvIG8EkabIiEoAkjaLtcjjoQXM4aEmazPsAJKlQJgBJKlQRXUCeBJakyWwBSFKhTACSVKgiuoB8HoAkTVZEAvB5AJI0WSMJICLOBX4feBy4HfiTzNzYr/quOL01/n/9PIC6LEkla+ocwJXAoZn5AuDHwAcbikOSitVICyAzv9VWvBY4sZ/1dT5Qod8DLEnSKBiGq4DeBlzRdBCSVJq+tQAi4ipgry6zlmfm16tllgObgS9Ms51lwDKAsbGxecUyqKFVJWmU9C0BZOYx082PiJOBVwMvz2meTJ+ZK4AV0HoofE+DlKSCNXUV0LHA+4GjMvPRQdXrkb8kbdXUOYBPArsCV0bEDRHxmYbikKRiNXUV0HObqNdzAJK01TBcBSRJakARQ0HUR/71E8FsCUiSLQBJKpYJQJIKVUQXUN31M1VZkkpURAJwOGhJmqyIBHDo2CJg65F/XZakknkOQJIKVUQLoL7c84RzVk0oS1LJikgA9XX/9bOAvQ9AkuwCkqRiFdEC8HkAkjSZLQBJKlQRLYCaR/6StJUtAEkqVFEJ4IRzVv3mUlBJKl1RCUCStFUR5wDqo/76PoC6fOlpSxuLSZKaZgtAkgpVRAugPtI/7syVE8qSVDJbAJJUqCJaAPUdwPXzALwjWJJsAUhSsRppAUTEGcDxwBZgPfDWzLy3X/U5FpAkTdZUC+DczHxBZr4Q+AbwkUFUuvauB1l7l88DliRoKAFk5kNtxacD2UQcklSyxk4CR8THgbcAvwJeNs1yy4BlAGNjY/Oqq77xqz4J7I1gktTHFkBEXBURa7u8jgfIzOWZuS/wBeDPptpOZq7IzPHMHF+8ePG8Ynns8c089vjmKcuSVKK+tQAy85hZLvpFYCXwl/2K5dCxRQDceOeDE8qSVLKmrgI6KDNvq4qvAX7Uz/rqq37qO4G9CkiSmjsHcFZEHEzrMtA7gVMGUalH/pK0VSMJIDNf10S9HvlL0lbeCSxJhTIBSFKhTACSVCgTgCQVygQgSYUyAUhSoUwAklSoyBydgTgjYgOtG8e2xZ7Az3sQTi8NY0xgXHMxjDGBcc3FMMYEvYlrv8ycNJjaSCWAXoiINZk53nQc7YYxJjCuuRjGmMC45mIYY4L+xmUXkCQVygQgSYUqMQGsaDqALoYxJjCuuRjGmMC45mIYY4I+xlXcOQBJUkuJLQBJEiYASSrWdpcAIuL1EXFTRGyJiCkvnYqIYyPi1ohYFxEfaJu+f0RcFxG3RcSXI2LnHsW1KCKurLZ7ZUTs0WWZl0XEDW2vf4+I11bzLoyIn7TNe+Gg4qqWe7Kt7svbpvd8f81yX70wIlZXn/WNEfGHbfN6uq+m+q60zV9Y/e3rqn2xpG3eB6vpt0bE0m2JYx5x/feIuLnaP9+OiP3a5nX9PAcQ01sjYkNb3e9om3dy9ZnfFhEn9yqmWcb1t20x/TgiNrbN69e+uiAi1kfE2inmR0T8zyrmGyPi8LZ5vdlXmbldvYDfBg4GrgHGp1hmB+B24ABgZ+AHwCHVvK8AJ1XvPwO8q0dxnQN8oHr/AeDsGZZfBDwIPK0qXwic2If9Nau4gIenmN7z/TWbmIDnAQdV7/cB7gN27/W+mu670rbMu4HPVO9PAr5cvT+kWn4hsH+1nR0GGNfL2r4/76rjmu7zHEBMbwU+OcX3/Y7q3z2q93sMKq6O5d8DXNDPfVVt9yXA4cDaKea/ErgCCOBI4Lpe76vtrgWQmbdk5q0zLHYEsC4z78jMx4EvAcdHRABHA5dUy30WeG2PQju+2t5st3sicEVmPtqj+qcy17h+o4/7a8aYMvPHWT1XOjPvBdYDk+507IGu35Vp4r0EeHm1b44HvpSZmzLzJ8C6ansDiSszv9P2/bkWeE6P6p53TNNYClyZmQ9m5i+BK4FjG4rrjcDFPap7Spn5XVoHeVM5HrgoW64Fdo+IvenhvtruEsAsPRu4u618TzXtmcDGzNzcMb0XnpWZ9wFU//7WDMufxOQv4cerpuDfRsTCAcf1lIhYExHX1t1S9G9/zWlfRcQRtI7sbm+b3Kt9NdV3pesy1b74Fa19M5t1+xlXu7fTOpqsdfs8BxXT66rP5pKI2HeO6/YzLqpusv2Bq9sm92NfzcZUcfdsXzX1UPhtEhFXAXt1mbU8M78+m010mZbTTN/muGa7jWo7ewP/EVjVNvmDwP20fuhWAO8HPjbAuMYy896IOAC4OiJ+CDzUZblZ7a8e76vPASdn5pZq8rz3Vbcqukzr/Bv78n2away3HRF/DIwDR7VNnvR5Zubt3dbvcUz/BFycmZsi4hRaLaejZ7luP+OqnQRckplPtk3rx76ajb5/r0YyAWTmMdu4iXuAfdvKzwHupTXg0u4RsWN1JFdP3+a4IuKBiNg7M++rfrTWT7OpNwCXZeYTbdu+r3q7KSL+AfiLQcZVdbOQmXdExDXAYcDXmOf+6kVMEbEbsBI4vWoi19ue977qYqrvSrdl7omIHYFn0Graz2bdfsZFRBxDK6kelZmb6ulTfJ7b+qM2Y0yZ+Yu24t8BZ7et+9KOda/ZxnhmHVebk4A/bZ/Qp301G1PF3bN9VWoX0L8BB0XrCpadaX3ol2frDMt3aPW/A5wMzKZFMRuXV9ubzXYn9UFWP4R1v/trga5XDvQjrojYo+5GiYg9gd8Bbu7j/ppNTDsDl9HqI/1qx7xe7quu35Vp4j0RuLraN5cDJ0XrKqH9gYOA/7cNscwprog4DDgfeE1mrm+b3vXzHFBMe7cVXwPcUr1fBbyiim0P4BVMbAH3Na4qtoNpnVRd3TatX/tqNi4H3lJdDXQk8Kvq4KZ3+6ofZ7ebfAF/QCtDbgIeAFZV0/cBvtm23CuBH9PK5Mvbph9A6z/pOuCrwMIexfVM4NvAbdW/i6rp48D/aVtuCfAzYEHH+lcDP6T1Y/Z5YJdBxQX8l6ruH1T/vr2f+2uWMf0x8ARwQ9vrhf3YV92+K7S6lF5TvX9K9bevq/bFAW3rLq/WuxU4rsff9Zniuqr6P1Dvn8tn+jwHENPfADdVdX8HeH7bum+r9uE64E8Gua+q8keBszrW6+e+upjW1WtP0PrNejtwCnBKNT+AT1Ux/5C2qxp7ta8cCkKSClVqF5AkFc8EIEmFMgFIUqFMAJJUKBOAJBXKBCDNU0TsG61RRxdV5T2q8n4zrSsNAxOANE+ZeTfwaeCsatJZwIrMvLO5qKTZ8z4AaRtExE7A9cAFwDuBw7I14qQ09EZyLCBpWGTmExFxKvDPwCv88dcosQtI2nbH0bql/9CmA5HmwgQgbYNoPW7yv9J6YtOfdwx2Jg01E4A0T9Voo58G3puZdwHnAv+j2aik2TMBSPP3TuCuzLyyKv9v4PkRcdQ060hDw6uAJKlQtgAkqVAmAEkqlAlAkgplApCkQpkAJKlQJgBJKpQJQJIK9f8BBaWMhtAGmnUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(X[0,:,0],Y[0,:,0],label=\"data\",marker=\"+\",color=\"steelblue\")\n",
    "plt.xlabel('X')\n",
    "plt.ylabel('Y')\n",
    "plt.title('Plot of data samples')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize the MINE model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "from model.discrete_minee import MINEE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "batch_size = 100       # batch size of data sample\n",
    "ref_batch_factor = 10  # batch size expansion factor for reference sample\n",
    "lr = 1e-4              # learning rate\n",
    "\n",
    "minee_list = []\n",
    "for i in range(rep):\n",
    "    minee_list.append(MINEE(torch.Tensor(X[i]),torch.Tensor(Y[i]),\n",
    "                            batch_size=batch_size,ref_batch_factor=ref_batch_factor,lr=lr))\n",
    "dXY_list = np.zeros((rep,0))\n",
    "dX_list = np.zeros((rep,0))\n",
    "dY_list = np.zeros((rep,0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load previous results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "load_available = True # set to False to prevent loading previous results\n",
    "if load_available and os.path.exists(chkpt_name):\n",
    "    checkpoint = torch.load(\n",
    "        chkpt_name, map_location='cuda' if torch.cuda.is_available() else 'cpu')\n",
    "    dXY_list = checkpoint['dXY_list']\n",
    "    minee_state_list = checkpoint['mine_state_list']\n",
    "    for i in range(rep):\n",
    "        minee_list[i].load_state_dict(minee_state_list[i])\n",
    "    print('Previous results loaded.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Continuously train the model. The following can be executed repeatedly and after loading previous results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd3gU1frA8e+7JQkQegKEjggqIoIg2Btg7+0n6lWvhY6KBUEgCSGgiCIdwXYVrwXLVazYGwoSpIP0FiAhpEAChGw5vz922bQN2YSETZb38zw8zJw5c+ad3ey7s2dmzogxBqWUUtWfJdgBKKWUqhia0JVSKkRoQldKqRChCV0ppUKEJnSllAoRmtCVUipEaEIPMSLys4g8HOQYThORZSKSLSKPBlA/XkTe8U63FJEcEbFWfqShS0TuEZFvgx2HOrE0oVdDIrJNRA57E1+qiLwpIpFlbKO1iBgRsVVCiMOAn40xtY0xU8uyojFmhzEm0hjjqoS4QpK/99IY819jzJWVtL2gHzQo/zShV183GGMigXOAc4FRQY6noFbAmmAHUVAlfXEpVaVoQq/mjDG7gK+BjkWXiYhFREaJyHYR2Ssib4tIXe/iX73/Z3mP9M8XkVNF5BcR2S8i+0Tkg5K2KyI3isgaEcnyHrGd4S3/EbgcmO5tt72fddt4t5MtIt8BUQWW+Y42ReQuEUkqsu5QEZnvnQ4XkRdFZIf3l8orIlLDu+wyEUkWkWdEJAV401s+TET2iMhuEXnYu61Ty9Dek97Xco+I/LtAXDVE5CXva71fRH4vsO55IvKH97VaISKXHeN1bSoiH4tImohsLdhlJSLdRSRJRA5445t0jPfyARH5vcC6RkQGishG7+s+VkTaisif3vbmiUiYt259EfnCG0Omd7q5d9k44OIC7+90b/npIvKdiGSIyHoRubPAtq8VkbXe7e4SkadK2n91nIwx+q+a/QO2Ab280y3wHA2P9c7/DDzsnX4Q2AScAkQCnwBzvctaAwawFWj3PWAkni/6COCiErbfHjgI9AbseLpYNgFhRWMoYf0/gUlAOHAJkA28UzQuoKZ3WbsC6y4B7vJOTwbmAw2A2sDnwHPeZZcBTmCCdzs1gKuBFOBMb9tzvds6tQztJXj3+VrgEFDfu3yGd7+bAVbgAu92mwHp3voW72uWDkT7eV0swFIgFgjzvm9bgKsKvG7/8k5HAucd4718APi9wLzx7lsd7/4fAX7wbqMusBa431u3IXCb9zWqDXwIfFqgrULvL1AL2An82/u+nQPsA870Lt8DXOydrg+cE+zPUKj+C3oA+q8cb5onoecAWcB2YCZQw7vM92HzfmAHFljvNMDh/dD5SwJvA3OA5qVsfzQwr8C8BdgFXFY0Bj/rtvQmxloFyt7FT0L3zr8DxHqn2+FJ8DUBwfOl0rZAO+cDW73TlwF5QESB5W/gTdDe+VO92zo1wPYOF3m99gLneff/MHC2n/19Bu+XaIGyBUeTZ5HyHsCOImUjgDe9078CY4CoInX8vZcPUDyhX1hgfinwTIH5l4DJJbxnnYHMAvOF3l/g/4DfiqwzG4jzTu8A+gF1gv3ZCfV/2uVSfd1sjKlnjGlljBlojDnsp05TPAn/qO14knnjEtochiex/eXtTnmwhHqF2jXGuPEcoTULIO6meJLDwSJxleRdoI93+m48R4qHgGg8iX2ptysjC/jGW35UmjEmt8i2dxaYLzgdSHvpxhhngflDeI6Uo/D8otnsJ/5WwB1H2/S2exEQU0LdpkXqPkv++/UQnl9H/4jIEhG53k8bx5JaYPqwn/lIABGpKSKzvd1HB/B8kdSTkq88agX0KBL3PUAT7/Lb8PxC2e7taju/jHGrAOmJotC2G8+H7aijR8ep+Em+xpgU4BEAEbkI+F5EfjXGbPLT7llHZ0RE8HT97Aogpj1AfRGpVSCpt8RzBOnPt0CUiHTGk9iHesv34UlCZxrPeQR/ira5B2heYL5FgelA2ivJPiAXaAusKLJsJ54j9EcCaGcnnl8E7fwtNMZsBPqIiAW4FfhIRBpS8mtXXk/i+TXXwxiT4n3tl+H5ssfP9nYCvxhjepcQ9xLgJhGxA4OBeRR+7VUF0SP00PYeMFQ8JyEjgfHAB96jzDTAjacPFQARuePoyS8gE88H19/lg/OA60Skp/dD+iSePtk/SgvIGLMdSALGiEiY94vjhmPUdwIfARPx9G1/5y13A68CL4tII2/8zUTkqmNsfh7wbxE5Q0Rq4umrPrqd8rRXcN03gEnek5pW74nJcDxdRjeIyFXe8gjvCdbmfpr6CzggnhO5Nbz1O4rIud547hWRaO/2srzruPDzXh6n2ni+3LJEpAEQV2R5apFtfQG0F5F/iYjd++9c7+scJp5r4usaYxzAAfz/TakKoAk9tL2B58Tfr8BWPEeRQwC83RbjgIXen8nn4bn8cbGI5OA5gfaYMWZr0UaNMeuBe4FpeI5Ob8BzGWVegHHdjae/OANPsni7lPrvAr2AD4t0eTyD52TsIm/XwPd4jiz9MsZ8DUwFfvKu96d30ZHytFfEU8AqPCdtM/CcjLUYY3YCN+HpOknDczT7NH4+e8Zz7f0NePqst+J5bV/Dc9ISPCd113jfnyl4Tg7nlvBeHo/JeE4i7wMW4el6KmgKcLv3Cpipxphs4ErgLjy/3lLIPxkN8C9gm/c17Y/nb0dVAjFGH3ChTk7iudRyNRBe5ItCqWpJj9DVSUVEbvF2A9THcxT5uSZzFSo0oauTTT88XR+b8fTlDghuOEpVHO1yUUqpEKFH6EopFSKCdh16VFSUad26dbA2r5RS1dLSpUv3GWOi/S0LWkJv3bo1SUlJpVdUSinlIyIl3lmtXS5KKRUiNKErpVSI0ISulFIhQhO6UkqFCE3oSikVIjShK6VUiNCErpRSIUITulLVyKz4gcwa+2jpFdVJSZ9YpFQ1kkojfTyEKpEeoSulTnqb//6FfTs3BjuM46YJXalKNDv2SWYmDDnh2507cyyzxgw64dutrubO/4lX5sw77nbmvPwO8978ogIiKh9N6Mpn3crF/PXbDxXa5s6tW5kd+wQzXxheoe1WF3sstUlz+h1HKSC/ff8ls+IHA5C08PsS670yehiveOsBbNsjpPofvwmA2YmP8sroYaVu//O3p5K6Z6dvfteOjSxc8GEgoZfoQGYmryY8zqGDB5k6Ko5Zo/z/bWxdv5IPpsXxyuhn2LLub1952t7dbN+6jtmxzxAfH88b454KaLt5uYf44OUROPLycDocvDI6jr+/fs+33Gl1FFtn5fcf4XJ4nqyYuXMHz40ez4oF3+B2uTiUmcGchBfZv2cXhw7m4na72b1/E2u3J+FyuMlIycFxpHD/2MGDB6nMIcuDNh56t27djA7OFXwvxw0iwhXJ1fc/zLtvz8NhcRI/pugzgctvZvwA9tKYOg4LtvB0Ius04MHH89uflvAoxuXk0TEzC623fMkf5ObmsvT7D0mTxjQyexlYpE5p3nnleQ7uSuPGfkOJad6cj96ZRcNGTbj8ylvKvB8rli4kpukpfDR7LHtpBEB8fDwA+1JT+WhmAm67nYGjJjM7diiIEzduUiW6UN1AbN64jm8/fJW8g3VwYeWA3UU0aYirDnut4b72ZsUPIpVomrCPFKIAGNT/EaKbNPNtr0WdI+RmWLj3yaeoW7++bxtHlzchDVy1uOLuPrQ/o0uhOGbHD2YPUUQ78xiUOJ7XnnuClEMNfYmvqXs/Me3as3P9TmzhOeCCdl0vIabVqWzbtI5ty5aRYq3pW9/X7uin2WOtRVMy2Y0npo4tarNhay4NZD8uE8agxEReHpnAfrsbgEZOBwMTxwHw3KhxHLEVTr42l51rrz+Xc3pcCcC2beuIrF2fqIZNfHXmxI5gtyWceg4LEdY8Uiw2LG4bbkv+A6s61xVuHur5+/zo5fGs3p9Hc4uFiDBhU27gJy/63HY/7338Fpd0u5Irrr+A7OxsNm/ezKeffsr1119Pt27dAm6rKBFZaozx24Am9JPcxJGJHLQ7iTZ7SZPCiSoQ02MHIAJOsXBR7zvpeuGlhZbPjB/IXhpRw2HjsN2JzWVj1NhRALwSN4gUb8K77LyurFz0DY/GzwBgwqhxHLY5sLvsOKwOop1ODILblgkWCxa3E3E24JZ+/fjo9fHc8fBI5s95CYfkMnjMLGbGebaLQGOzjwFjpvv2K8rpYnDiWN6ZOZ5zLupNh07nAvDR21PJzckhe3c2qbZw4uPjSd66gfn/mcxeaUSE005ugURytL3ZsUPZY/E8x7mRLYu9znrFXqczWjUkOz0NZ3Y6KRJNtNOB1bofgP5jppOelsaShd9x9c13M3Z0Ii5r4afixbiy2WOt7Zu3uK24Lf4TjNVlK7Z+Y/bR8/8e4tv/fkIOEYX242j7/ca+5NmfuEfBHcYea6RveVMy2E0Dv9sLVA2nnTDjIsdiwWV1Eua0k2crflQM0NDhJt1euAOhpe0AO5x1jrmNSIcVu7jJtHnyWrjTTuNWhvDMbDYeqlFqjLWdVm68sQer/vyblZm5Ae7ZscXHxzN+/Hjy8vKfoT548GCioqLK1Z4mdFWi+LgxIIYmJp0UaQjAKY3t3DdgpN/6z40ah0tg1FjP8qLJv4nZR46zMRFyBGNNxyrCXhphd9pxeD+8R9cpuO7RI6UYk84RcZHhPQoukRGQwP52Y0wGe6RwMopx7/clYYAo0shwNymWJOs5hCx7ydsRtxVTQmINVLjTXuyI80Rr4DQckDC/3Q6q4nXu3Jmbb765XOtqQleFzI59wpsQ3YWS2lHittJYUjjijCbTXri8YPJqYtJ8R9gliXFnssdSv0hZFm53DVJt4ce3I0pVU2effTa33FL2rj84dkLX69BD1MxRz2IRB/3HTmRW3GAsWHCIg8Hxs9hjOfbPVmNxkUI02IuXF1RaMgeKJXNPWT09Ha9Oana7vfRK5aAJPUTttYUBYQCkSn5f3bS4ASCNgxSVUgrAarVWSrua0EPItPFP4so9jOAGi+fs/oxRo8CW/zanazJXKug0oatS7T/coNhJrTSbvsVKVTUiUintBtSTKSJXi8h6EdkkIsXuAhCRl0VkufffBhHJqvhQ1bFMjeuvVygoVU1ERkaWXqkcSj18ExErMAPoDSQDS0RkvjFm7dE6xpihBeoPAboUa0hVildfiiUtM5w8W5PSKyulqoSzzjqrUtoN5Ai9O7DJGLPFGJMHvA/cdIz6fYD3jrFcVSDngfQSb85QSlVNlXWVSyAJvRmws8B8sresGBFpBbQBfixheV8RSRKRpLS0tLLGqoqYNjLWd2u5qr6ipfyfhdqOyjm5Fogazvyk1NS9P3/amsVtN15L24YFlpOFxR28WCtDlKlZ4rLzOl3BkCHFB2Vr1Mhzw1x4eOXcgxFIQvfXe1/S3Uh3AR8ZY/zeOmeMmWOM6WaM6RYdrYnoeMyOfazYrdHBYnHbiDHpJ3ajRrC6bNhcnqRhd9pp6s6gjuPYr0mkw0aM+wBRztLv7mxW89i3flvcNhqbwJJxE2cuLepAHYelWGIbFDeD+Ph46jn8nygrmrQL7uNZ553im27qzvS7flN3hm86xpXtm4525tHIdYQY936aujO4v89thdZr7Dz2/j+TOJKm1gwaOA33PJ1AjDubFhEH6Tt6Mmed051/DRlJE9chmkoWfeMnE5sw+pjtgefu4ZYROYXmS3P0b6CkeYAI75ePxWWjhRw+ZnsD/nUrI0c8w79v782/rr3IV97Sahj+zEDEbSHSFUHf4fmDoT3+yEM88/Rw7r79AeLj47n61kuoUcMz1EDHjh199fr168eoUaNK3afyCuQSiGSgRYH55sDuEureBeiYnSeAvxt2TqSC45o0YBf9xswCPB/Ahk43OYT7bmc/eodpmNNOA9lPitV7ZGMsIG4amFTEGU2G1UZ9SwogiC2C+g2bkJG6jQw8X/6NzT6sxkqeHGHwmFcAmDHuSdIctaljjtA3YaovBgzc/687+fbTd3BkWX3jl0Tad9Mv3jPI1+zRT/rGRyl61+vFXdvR84Z7fAklWtKxuiHP1YAMm9DQ6WJIomfZjFEjSbPZwVhoVuMgeXm5pLkLvz/9E58vNP/rgs9JTdmGw5E/vkejljXI2nPIs6/OXFJtEdRw2Hly3EhmjBqFzZpF1ytvptuFvXxxXXnDPaxcNJYIyaNv4pTCCdAIDV0u+iZOJT4+ntoOK/3GvcTzo8aRa3Ngs+TSL2FSobg6NP2SnKx91I1qwm0PxjN+1DjybI5i48Mc3U7f0VN9Zf0SXqKo/mNfKFZW0rANVpcnHT04/MVC+1E0qcfHxyPGQqdTa7FiczZNwjLJyK3HIbsnvjBjKDySDfzr4Zv5bu5b3D08nrCImsXajHZBjx5tcTkcNG7bCYBWHS/0LPzqd8Rt4cH4WADiEmKLxV6vmSdFtu/Y2ldWs2ZNBgwYQIMGDVi9erVnHyvpcsWjAknoS4B2ItIG2IUnad9dtJKInAbUB/6s0AhV0FjcVurIbrKKXLvexHWQQ6YOR4/f2p59gW+Zvw9iXMJoZsUNplG7s7jt/pEkjk7EaXXSkD0ADPF+GZRkTuyjuE0E/cdOL7Zs0MiXmBI3kC6XXusra2j2Ymx22pzagX5PjS+2zlFHB6M6auLIRASobcmi5w33ANC81hFsYeE88Ni0EtuxWDOBRtRwWXlk+AQAvv7oLbZvWIE5kovFhBVb55KrbihWdne/YcyMH4jVXYN+iS8xO+5RGp/pOXk2KDGxUN0aDjs15QgAT40rcORrhJpOK8PGFT4K/Ne9dxLVMAaAy6/uwrIFv3D1vx8pFsOdfQtfxNakdibZByK5a/D91KpVlxcnTSq2Tlk0tWZRt0kzcg/lsCstjDybg7oO4YjYiKqR/yuvhWU/Lqf/9HRKnSM0a38mPa+/l6M3zyfv2szCBR+wbkce9ewHuPmuO/np3Q9JMbWo44Jmrc/kgdH5Xy4daxlWH/T8ImptcfJAfKKfLXkMvO927DX8D+zV7847OLh/v99lAI0b53926tQ59h3aFSGgsVxE5FpgMmAF3jDGjBORBCDJGDPfWyceiDDGBDTwtY7lUn7TR41mn6383/SNTCp7A7jBqKFJZciYWUwbOxTjyiOqeWv2bd/CTQ88xadvTyRTmmB12Rg91v9PyNmjnwLLEfqNKZwME0ePw2l10IBUHo0/djKvDnwjSjrtPJPof1CzUJBz4ADTX5hGfete+o2ZctztzUl4nN3uejR25jKgyC+Y8tqwPonWrc8iLIA+6h/fmswfm3Lo9/DtRLc6vUK2X5LMzEwiIiJ83TDHQwfnCiGTYweQZQn8bs8G7hQAMiz5lzXWN2lkFuheiDEZ5OJCcGOx1cR6JIy2PTpx5Q19Smx3SvxAMmmExWUjtoSEXuI+jBxDlt3Q0JrFkNGTy7RuVbRh7Sq+fXceEnaAQbHHn+hOJnMSHufWgc8SFVXK6JrKRwfnChHff/lJmZI5wKMJnr7mo90f9cjgsTEz/A5hWxZRTVqTmXKIKHOozOvGnBKFbftqBofA0TlA+w5n0T6xcq4rDnV9Y6v/F3pVogm9Gvnrz3UBvWMxrhxyJQe3WAuVOU04gxKn+uZdfvp2A3VP/2HMf+9VOvfoWeZ1/+8hPW+uVGXQLpdqJNCj6RiTXqzfWikVGo7V5VI1LmRWFcLutNPYpGkyV+okpV0uVdjk0f3JpjlRlj3HvCM0yuni2n/dyymnnXYCo1NKVTWa0KuwCGqTZXWSiv9k3sBhcNn2MjgxNE4uKqWOjyb0aswueTxayk05SqmTh/ahV2klD4Lf2HmElmefeQJjUUpVdZrQqzBjSr4b1GHJ4ro77j2B0SilqjpN6FXYQVMr2CEopaoR7UOvombHPkmOvXax8roOCxFykCvvfSwIUSmlqjJN6FXQbz98wx5L8WQOUMO2l/5jio86qJRS2uVSBS395bMSlzml9AczKKVOTnqEXgWVNBhDE/bSP0QGtFJKVTxN6FWRKwYs7kJFN1/fi87dLiphBaWU0i6XKmdybH/22wsnc4vbpslcKVUqTehVTIQU/9HUyFLSI1yVUipfQAldRK4WkfUisklE/D5iTkTuFJG1IrJGRN6t2DBPHmKKvyX9vQ81VkqpYym1D11ErMAMoDeQDCwRkfnGmLUF6rQDRgAXGmMyRUSfJ1VOeywNCs03cR0MUiRKqeomkCP07sAmY8wWY0we8D5wU5E6jwAzjDGZAMaYvRUb5smr/9iJwQ5BKVVNBJLQmwE7C8wne8sKag+0F5GFIrJIRK7215CI9BWRJBFJSktLK1/EIWxy7MBC8w2des25UipwgSR0f0P+Fb1U2ga0Ay4D+gCviUi9YisZM8cY080Y0y06uuQHNpyMpiY+RZalcE9VuGV/kKJRSlVHgST0ZKBFgfnmQNHLLpKBz4wxDmPMVmA9ngSvAmTLLf7A5r4JU4MQiVKqugokoS8B2olIGxEJA+4C5hep8ylwOYCIROHpgtlSkYGGOqexF5qPcevRuVKqbEpN6MYYJzAYWACsA+YZY9aISIKI3OittgBIF5G1wE/A08aY9MoKOhRl2Av3bPVLeDlIkSilqquAbv03xnwFfFWkLLbAtAGe8P5TxyncaS+9klJKFaF3ilZBkTWzgh2CUqoa0oReBQ15Vq89V0qVnSb0KmByvD59SCl1/DShVwHizL/dX9z6liilykezR5AtX/I7mbb8+7Qai46aoJQqH03oQbb4iw+KlPi7MVcppUqnCT3YitxQdNjo+C1KqfLRhB5kxlgLzQ9N0LHPlVLlowk9yHLdkcEOQSkVIjShB1mWPf+EaF3JDGIkSqnqThN6FTI0bkqwQ1BKVWOa0INo+qhRwQ5BKRVCNKEH0T5b/thoTXW4XKXUcdKEXkX01eFylVLHSRO6UkqFCE3oSikVIjShB8nsifHBDkEpFWICSugicrWIrBeRTSIy3M/yB0QkTUSWe/89XPGhhhZnTv4gXM3rFn9AtFJKlVWpj6ATESswA+gNJANLRGS+MWZtkaofGGMGV0KMISlNGvmmHx76bBAjUUqFikCO0LsDm4wxW4wxecD7wE2VG5ZSSqmyCiShNwN2FphP9pYVdZuIrBSRj0Skhb+GRKSviCSJSFJaWlo5wg0Nk2IHBDsEpVQICiSh+xug2xSZ/xxobYzpBHwPvOWvIWPMHGNMN2NMt+jo6LJFGkJqSX5Pl81lP0ZNpZQKXCAJPRkoeMTdHNhdsIIxJt0Yc8Q7+yrQtWLCC1EFxkBvpE8oUkpVkEAS+hKgnYi0EZEw4C5gfsEKIhJTYPZGYF3FhRh69ljq+KYPm7wgRqKUCiWlXuVijHGKyGBgAWAF3jDGrBGRBCDJGDMfeFREbgScQAbwQCXGHFIeGzsr2CEopUJEqQkdwBjzFfBVkbLYAtMjgBEVG5pSSqmy0DtFlVIqRGhCP8E2rlnjm67r8HcBkVJKlY8m9BPsf+/9zzc9dFxcECNRSoUaTegn2CGbM9ghKKVClCZ0pZQKEZrQT6DvPn3fN93YeeQYNZVSquw0oZ9AG5b94pt22bKCGIlSKhRpQj+h8q9qcaJjuCilKpYm9BPI5cofkOzx+KlBjEQpFYo0oZ9AGTa97lwpVXk0oZ8gBW8oqu8IYiBKqZClCf0E+eq96b5pkaLDySul1PHThH6ChFvzX2qXnLxPa1JKVR5N6CdIiuSfEB06dmYQI1FKhSpN6EopFSI0oZ9gdqdef66Uqhya0E+whqJ3iCqlKkdACV1ErhaR9SKySUSGH6Pe7SJiRKRbxYVY/U0ameCb7j92YhAjUUqFslITuohYgRnANUAHoI+IdPBTrzbwKLC4ooOs7g7Y3cEOQSl1EgjkCL07sMkYs8UYkwe8D9zkp95Y4AUgtwLjU0opFaBAEnozYGeB+WRvmY+IdAFaGGO+qMDYlFJKlUEgCd3fACS+Wx1FxAK8DDxZakMifUUkSUSS0tJOjptrCt7y34i9QYxEKRXqAknoyUCLAvPNgd0F5msDHYGfRWQbcB4w39+JUWPMHGNMN2NMt+jo6KKLQ9KP817xTVtrRQUxEqVUqAskoS8B2olIGxEJA+4C5h9daIzZb4yJMsa0Nsa0BhYBNxpjkiol4mpmjzT0Tfd7OuEYNZVS6viUmtCNMU5gMLAAWAfMM8asEZEEEbmxsgNUSikVGFsglYwxXwFfFSmLLaHuZccfllJKqbLSO0Ur0cy4Ab7pJiY9iJEopU4GmtAr0V5p7JvucNG1QYxEKXUy0IR+glzS+5pgh6CUCnGa0JVSKkRoQldKqRChCb2SzIkd6puu70wJYiRKqZOFJvRKsttS1zf9WOIrx6iplFIVQxO6UkqFCE3oleDT918PdghKqZOQJvRKsPyf/NGGG5rUIEailDqZaEKvYL98+2Wh+SFjZgUpEqXUyUYTegVb8+vCYIeglDpJaUKvQG9OH8deW5hv3u60BzEapdTJRhN6BcrbW3gArqi6+nhVpdSJowm9Au2x1Ck03+/psUGKRCl1MtKErpRSIUITegWZFDug0Hzn01uUUFMppSpHQAldRK4WkfUisklEhvtZ3l9EVonIchH5XUQ6VHyoVdsBS/7Y5xi4+a6HgheMUuqkVGpCFxErMAO4BugA9PGTsN81xpxljOkMvABMqvBIq5Fwl17dopQ68QI5Qu8ObDLGbDHG5AHvAzcVrGCMOVBgthZgKi7E6qeBZW+wQ1BKnYQCeUh0M2BngflkoEfRSiIyCHgCCAOuqJDoqoFv/vcui1ZsKFTWL2FKkKJRSp3MAjlCFz9lxY7AjTEzjDFtgWeAUX4bEukrIkkikpSWlla2SKuoLUtXBzsEpcrk7dlvsvDn34IdRpWXvm8/K1dsDnYYZRJIQk8GCl6y0RzYfYz67wM3+1tgjJljjOlmjOkWHR0deJRV1KTYAYXuDAVoYvYFKRpVXu+9+TrjYieUWm/imLHMnDiBpN9PXDJMjJvEgvkf+ebnzprBz999QWZ6Ore+/TPx8dNLXHd/Zib7UvIfrjJhzIs8H/8Sw9p3YWR66T/OP3n3fUYmvk6Tn5YzZULgvzpnTbuRn4oAACAASURBVCl9/P83X/sPkybODKi9377+hj3J233zX877jDdmzA04Hn/27dvPBZ8uZtqrX5ZY5+4fV3NlRvYx23HkOXDkOY4rlooUSJfLEqCdiLQBdgF3AXcXrCAi7YwxG72z1wEbOQkUurLFq/+Ykj9gVd1zo0bjsNQhNuHpgOrf/fr3bGhUl6Qbzi3X9p6Y+AF1s3cTlzDU7/LnRo/lmht70fnc831lu3fu5O3Z/6F2LWHQiPwfgk1+Wg5A31//RDCsbXU6p+zZxfcdT+fBtd+x2XIKFzQ+xO33PQhAQuzL7Ipuy5hbu/K+/VQWX96VmBnTuPaW2/jsg/fpN/SJQrFMiB/Ly5fe4JlxwL0T3qVT3UyatmpDVloKn+U2o7d7E/f1H0Ty9i28+fonbGrWkQXtm3Drym3MfOxmEmInsaXZGdxQJ43LenkeGt4wOprHXvqQD85px4jFC3hs+DOk7NrB9198wYHsg0y/rDcAt075lGYZm5l2WU9PDCt3Qot6/NHiIuK9Mb4+eRqZ2YdJrtGa97u298X+9Jx3+PDsi9l2SS9f2T8NatDkp+XctHoHX53RghouQ9ddGfzUJopuew7QO3kZH7Q/jy0XhgOwPqKlb91nx71J+/pOzjjzNKZvOEyzfalc17UpS5auIdMVyasXnseO594ixn6QVTVaMuSCGJKWrCAysgZ33t0HgBFtOwNgJszi12Yd6LpjI/fcdj5tTzuTw4cOMfHl//LEkDuJrFOXOyKawMZMvt6xiy4XXMBD0a0gGp71vucPJK3joau78vbXywA4YrFwSYvaAFzQszvp+3KYtmA54x/qzemL/6FLSjYDG4SxpW4Eb1kb8cjhI4x59Rv63nAubdo0pc9bP9EoN48V7fM/39u2prBmzVYWZhzi/06L4ewenutCLv5mGdsj7Vx02MLjp8Zw4WmNAEjJOMTibencdE7+sbAxhs/2ZtGjXi2ahNkR8df5cXzEmNLPX4rItcBkwAq8YYwZJyIJQJIxZr6ITAF6AQ4gExhsjFlzrDa7detmkpKSjnsHgmXZoh/57JtfC5XFx8cHJ5gC1vz9J59/8iWZNU8D4L5r23Nbio2scCvPr/mRS6+8npVLlrBq3Q5GJY7wrTdi3H948wLPhyzlcs//F3y6mDP2pNM4I5WsyPosOKMlbfbnctnKX7C6c5nS6xYAHvhzBb90OJ1OO1M5Petvho6OZ0zsy7hsNanp2E2nLmewZGUq53ZuzMJVOfx12hlctup3pl1xld99aJbj4MqVa3jzgs7Ucri5cn0y/+vY0m/dfr/+weo2Z7CwRf2AXp/+v/zGZ93OZ08tz7HMVRtSWNC+SbF6XVOyMUCPf/5mf2Rj3u12eqltR+W62Bdh9bvM6ja4LBX/AT7qoYVLef3CrpXWfjA0Pehgd63gXzFmdxkc1vK9d2eJjVXGCUDzCDvJuZ6j+ZGnxDCkVfEDwkCIyFJjTDe/ywJJ6JWhOif09atW8d7HHxcrr+iEvmLp77zz3XYubGu4+Y57S6w37bkx5Bx0MiJxLLe+/TN/tKgXUPstsh3839+fsqz5pfzQtlGhZY/+8BlTe95UwppKqeNxR5P6TDujVbnWPVZCD6TLRRXhL5lXhv9+t425PTryX2NYNzIOl6UWH/e4gvQIK3neI4Zwl+HIeZ7EO+Wn5RBgMgfYWdvOi5fe4XeZJnOlKk9tq/9fcsdLE3oZzXx+BBBeuNBAjXLeTHTfnAW02JtMjdwUrOYQU3oVSLA9OgLgFvF1bxR1pJw/BZVSwWOvhP5z0IReJi+PTGC/PbxYefyY+FLX/eX7r/jrl994euxzDJzyKZ90ak3XlGyWtmsM7crXl6aUqp4ublC7UtrVhB6AqSPjyLALlPP8zEuxzzLx8jvhiv/jYPw0Prn0YgCWNqmcN1UpVbX1alin9ErloKMtHsM/K1cya9QITzIvgc1PV8v4kaNp8tNy37+Jl9/pW/aKN5mrklnKeKL+vOT99PvldwBuXLOzlNpQy+EutU6nvQeLldU/4uL0jMP0SVpPzEGnr/ya9XsAeHDh3zz209c8/esXPLBoJQANcl303Jx/E13HtPx2u6TmlBpHWbQ6kOebfuGfxaRc3tl3xVJZXb82+ZjLn1j4M2+kb+LOZZsKlZ/tfd0u2JlVqLz/wr9IubwzQxcXfkTjfYvX8K+/1gJwesbhErf30OLVRB92Fiqbsm29bzrmoJPLtnoeMBPu8v/3c/OaXQxavrVYub+/tyc2pXBGkXiu2J5B8kUdfe/953Yry9o0Zsu57Xmjfv59NTcePnb/+C2NAj/PVVZ6lYsfs+IGke5uitNa+g0DTWrl0f/p8b75hLiX+Lj7paTWrJ4/fh75dRGvXnIeALes2k6DA/uKXQ7X6JCTvQX2r2V2Hv9e/yO7D9fn1YsLjwox4OefMWIlN7wu/zmvU6FlV25MpVnaHuod2cqWRl1puD+d+nlb6Xn91Xz3+QK6XXQuPa+8nsTYF+lzz1W0Pe0s3/XmI//8goy8Osy69BL6/raIhNj+hdp+8oX3OBxegxb7VuOQGmyLOYMlrRozaN33nN6hEwv/WIFLwoktcg38o5M+Zl6Xtryds5MlS9ZyxFqH2iada2+9njM7FU6OibEv8k/LTozv1ZaWrdse83Ud8vInNM1az4gxIwqVP/PcXN467yz6JK1nf2Rt/olpyCXr1vP8iPsYHzsRxLA96gw65qwj1Vmf1y46l/6//E58/GCefv6/zO1xJv/+czlRrr08OXoYPT77i/M2bWfKk/nnYkaM/w/1HfsYFvcUAJ/89wN6XncVmzds5pzuXXn3jbnUqlWTFf+k8MdpZ9Nz0xKeHpX/umzdsJEjjjxe/3w5c3ucyT1L1vHSsD6F9iMnOxunI496DRoWKp8+aRZH8tw8OXyQr2zuW3P56VBD4q47k1YtC1/pcfT9TVyzBBs2hp/ZBaDQF9Ov33xPVKOGdDinC2/MmMtV115Mszatef8/n/F4q1Zc989uvjy9KQD3LNtMuMvF+KfzD6rOnb+EnbXtPPnPDvrefQV160YCsHTJP/y4YhsvtW3C2043V/Y+B4DUlHQaNym8X/68+vtmrjqtMU0b1CQ3z8XSw7n834rN2EV466w2XNagNpYK6DvXyxYD8ErcIJp16sHWFX+SSQzGUvpRXIw7g34JU1m6+Ff+9+16rO48Zl9yYaXGOeCnH5l1ef5QOT127Wdxs7rF6nVOzeGl0+2c2bkHU8fF8Vn7q7h09Z/cd++NTPhmHS1S/6Zd26as2Z5Hg7BMXj/7GlJr2vwe0R39kC1uWZsVfy/m7HN6MPu9X/j71PYsbxzJU798wlPxCQAkxr7A/prNmDj8nmLtvJwYh8sJO2t35INz2jFuzQ88NPjJMu3/bz98w6qkvxn4zLMATE5M4PFRsWVqoyrZvWs70//zI/++owft2pc+6vSsF15kwDBPYs7KyODrTz6mz8OPVHaYPuPGTuPRx+6ndp3K6TKY+Pws3MAzwz3PF7jwf4vZXC884F8ar0ybx519evHxBz8QXb8WN999bbE6Sxb/w8rV23jooav9tpGRcYAGDSpn/yqCJvRSTBs1mnRb4JcRRTldDE70PF5u/apVXLrPVeEx3bpyG590ak2E07CtdxcS4l4iK7I5k57+PwA+fvtVtm7exlNjxhVa7+Wxozl8xMqzifFl2t5/X5vDoZz9PPJ48btE48bMxGkJZ9zowmO879iyhZSUXXS/QLuRVOVI+uU31qzezP2DHgh2KFWGJvQCJsUO5IClEXannbqWZPb5uX2/NH1uu41vv/qUyV2u46D9+E9DnLv7AJdu/oa8vAgOhcdwRpM87nlkEBNGPovNbuHJ+MTj3oZSKjSctAl97bJlfP/pK/S6ZQAffvI5xnL8R9INHIYsayQzL+9d5nUH/PwzLks4W5q2peP2n9na+Dx+at+cDVedc9xxKaVODifdnaKz44awRxpic9lxWmOY97/PwHJ8X1yRDitPjRvNxLgRzLys7Ml88A9fMipxZIGSK48rHqWUKiokE/oe8ZyR9l2lIuVP5lHsxSU2Hhs3ldiE2cy57P+OWf+M9MMMydvAkbzD5OTksCdNqBt2gMcSx5Q7BqWUCkTIJfRXXhjFce+WEZq4c8BymP5j8sdsnnNxsQc1+Vy3bje1D+Uw8OoOtO94//FtXymlyiHkEnr2/giwO0uvWIIYVzb9xr5UrPz2t36ClsWHaL1u3W7a7v6rzFeVKKVURQuphL749+84WMZkbnXZqF8rB8mx4bbuo99Y/09b+d1PMgd4feC1QPFrXZVS6kQLqYS+8OtFAY+3Ut8pREg2hznE4OHHfhTWUxPehe6Fb/ro9+sf2F0HoJy3ViulVEULqYR+wF74ssT6DjgsdnJtDjAWmpCGGxdHjJvHEkt/7uFR7xRJ5l1TshkTN7BCYlZKqYoSUgm9oCbspf84z5H3tLgB5Lmh/9hZZW5n+Pi34PyzffM1nW6+7KN3Riqlqp6AErqIXA1MwfNM0deMMc8XWf4E8DDgBNKAB40x24s1dILUcNrpn5jfjTJkTNkT+VH/KZDMAR76ZQH01huBlFJVT6n3rYuIFZgBXAN0APqISNFRhJYB3YwxnYCPgBcqOtCyuOzq7pXS7pAfPmNk4ojSKyqlVBAEMhBJd2CTMWaLMSYPeB8o9MBJY8xPxphD3tlFQPOKDbN0Mwr0afe4qOx3cpamfWYuIxPjKrxdpZSqKIEk9GZAwacGJHvLSvIQ8LW/BSLSV0SSRCQpLS3NX5Vys1TCszpGjX3VN/3rredVePtKKVWRAsmC/kZk93svvYjcC3QDJvpbboyZY4zpZozpFh0d7a9KuaXiaa+Oo2xP0/7h60+ZNW2s32WvXXTuccdVFU2YcjPPjx0c7DCqjTWr/+bvv/8E4NdfvuTtt7sy/rm+heo8//LtjH/pTn+rH5eU5B1MeWlYhbdb1PzP57Loz+8rfTuqcgVyUjQZaFFgvjmwu2glEekFjAQuNcYcqZjwysA7XkvRSxePen7yrThzaxPmqEf95ovJTO7BsNFT2JsRT/szs3nl9W/IymnE8Mf+B8CEUcOh512+9Z+bdAe39BrN6Z06+W0fYNyYh2nQ7B/Sd3dgVOwcZkwaSXZ2FsPjZgAwd85L7Ni9AamdSd06KTjdNizixpnejaEjXvS189qbPWjTah9hZga/Ln2HWu4Y3HVW0KllX3pdfzsA774+mYZNZ7BqczvyjkTiPliXMGstho2cTGLiw5zStAsO25scyI1kSN+ffW2PTXyQCy5YRcbezb6yxIn3YA87iCM7ioh6u3lq8FdMmnkVR9LbMmJ0ydfoz//kLdbu/B8WDLVqZhIVmUVU44P8+cdljBr1uud1GzuA7hd+y99rzgAg3H6IR/v/yMa1y/j0h9F0Pe0hrrjyFl+bkyc8RYuW7dmW+Qm52VGMfOYdAF4c/zgO5yFGxM7x1f35h69ZmPQuZ566lA3rezHs2amF4nvhuUcZNsJTNmH8YJy2LMJrpWMRN1HhV5K6bwfdu/YEYPXqxQx6vPDY8kdt2XgPNevm8fz4XkQ3XUbrlllEhC0qVOfcs5cB8MqMOJo1+4ADB2twzz3LmDjB88UZ1WQxKbs7073HjyQn1+f++zwjje5J3sa3P95GckoHRg6bC8APCz7i8KFcrr/lXuZ9eT8du+xg/IR92CSSpo3a0Si6GVded2uh7X/4/mzuuKsffy78hfMvvLTYPrzw4iCMcfDM03N4b+5kklM3QPgOnh7yBQC1asWzd08YsI4PP5xFZmYqffvGF2tn3T9LWbdxOT3O7c3aFb9x/gU3Elm7NhMm3UeD+hvIyjqLc8+6mmbNT+WLBbOxip1HH53CoYMHWbJ4AZdecWuxNmfPGc1BZxJRNS/jjjuH8NbcRGqG1+W+B55h1msX0apxGpF1Z/DXqqm0ariVPEs/7rnT87pOf20YtcN+oFWbsWQeSOWW6/7t9z08FofDwXtfPUefa0dgt9uZ9cHdXNBjKNtSVhJ5aDwp5hEOpPxO61a9Oa/TndSvE+Nb9/15N1ODvYRHZHPldcvITNuAPTySOvVbkpm2ifrRp5KXl07OwQ00qH9+mWMrq1KHzxURG7AB6AnsApYAdxtj1hSo0wXPydCrjTEbA9lwRQ+fGx8fD0AjVxqNOv5K/ehDpO6pTYP6h0hacinnX/xjQO3k5tjYl1WbRa6+/K+15+TqUDOBbvyF47AVew0XW7ZFc0prT5fRwf1h1Kqbx1/Lu9K981JfO5u3NaJt670Bx5+VHkG9hrnHrLNmY0sOHa5Ps4bbaNpsv986u3fVLbZs8bKu9OiylH82N+P0trtKjcXtFCw2U2Cb9ahZI4sz2+3w1UlO9tw527x5pt82Vm1og8NRg3POXFvq9ooybpACvx03bImh/Sl7fPNHDtoIr+X/juClqzuQlxPFeT1+RQQOH7CzamtHunsT7lHZmeHUrl/4uCNldx1sVhdRjQ/iyrPw17ILOb/Hb2WOv7wOZIRTp0F+TM48C7Yw/0/OWrb0Jrp0/czvslXrT+HxAd/xnzcm0qL1se+3+GvphVgjMunqfZ/+Wt6D7p0X+6179G+9qCUrunHu2YF9lpeuOpeuZy0BYGdyfTJyGnL26ZtKWcu/9L21aNio+LNfd+5sQIsWGYXKNm9vTNtWqX7bSdlThyYxB8q27awLaFjvjzKt07LFQ+zY+TpWayQXX7QYqzWiTOsfddzjoYvItcBkPJctvmGMGSciCUCSMWa+iHwPnAUc/dTtMMbceKw2KyehGy6+5J0Kae8B3sMhYQD819xWIW0qpRTAqW2H0apVv3Kte9zjoRtjvgK+KlIWW2C6V7kiqyB//vYtQIUlc8CXzJVSqqJt2vxCuRP6sVT8pSFBsGbZkkprO8FU/gkppZSqCCFx639uxi7CI2pUWHvrOc03HU3FXl6plFKnnPJEpbQbEkfodncYdeqULfH+vcZzs+vh7OLDMybIeN90HYqfLHHk5l8aeSAjvNCyg1nFu2qK1gnUot96szPZ/7C9pdmXWotFiy9i7aYWbNkWzfrNTUtdZ/2WGLZtjypUtjclslzbP9Gy9hX/Qt+zbQiLllxAjp/3pDSbtpb94eHgOSlbEZasPLv0SgUc/XsGz4n90mSm1STXz99+eRy9gqmoNRtbVUj75ZGVnn/CcdO20v/2S+PKEw5mhZGbYyM1td1xtdXu1JG0bjXguGPyJySO0I24aN8+/4xzzys2kzjmESIa7uapwV+SGN+PCFskT43yPLhi7fLlPD3kGMPe/rQcgHv/WkvPZzYz7rn7MU4b2JyEu+rz9KjJTJ9zKWecmozdNoqklfMQixunI4Jnn5zna2bSzKsItx9m0CO/8vzLt2IPP8SRw3VoaD8Duz2CdOcfHMlogrV2FnVrp+A2Fg7l1qFNveu5475+9LwiP6RXXr+Qdm1S+POvixg1/K1C4U5MfBxHzd20bryBxjHZLFp8ISNHvB3Qazd7yjgO2n8hMiKbmDoPcsudD7Jt0yZSUrZz3kU9C9V945Xniaz3Lg0bHWTRkvOpF34Kgx5P4K8/f+bHv6bQtMF2mnmvsFm3qTkHDsRgteXSrdMqwHPly1nttwKwfUdDLu4+m3kfz8JqsTJ8ZPHxdma8egmnt91F0sqz6NL2XpI2zMNmz2XYo/P97suLM67BkdUKaqRz7WVP0/NBz1VKn30yly3bZuJw2TiwtxP3P/Asy/7+g107/2G/WYgjL5LTmq9jz66bOXQkhzzLXkY+/Z7n9ZkWx77s3Rhx0CyqA1n7Uxn6lOfvaOK0G3DmNKdeZEPCw2vwYN+Rhd4zf8Y/9yB1Gm5mcN9fPDFPGILTnYs7bD84a3BKTDd2pKxn+LDpvP/2NLak/c4p0RfRuevF/LnoB/bu24IbB93P/YFD+8NYubkLtW0tefrxQsMrMWf2WGrVqEvXcy7hm9+e5qzTtpCeMgy7Xcjcn8aDD48sVH/S5MfIc+Uw/EnP5aZTZ19Os/pp3Hbnar75qj32CBdLll/CuZ1/JTfHxnU3rvet2/MKyMnOxunIpV6DaGbPepZbbx9MzyuasiLpZ77+7VXaNVtBvajDvnXWbGxJniOSLh3WetvwXEb76quxXND9Bn5Z/AXhYV9zOK8mFnGRm9eFJwZN5aWZDwH7qBW+jzM7xJF32JMY81zPkZ2zj1qRDeh1yS2Eh+cfRPUE8vLyWL7sRzZsWUjO4T107nY3K9Z9xk09h/HxF/2JCt/Ovrz29L37Pez2/C+61TuWkZK+lSu7FL7cMiVtAz999yC9r5rLihWeU4g1XbdTr8GZOI4coNMFnssqMzIWkp2zllYtHwEgJ2c9kZGnUWmMMUH517VrV1NR4uLizPc/nGK+/+EUM2lmr+Nur/GPy0zjH5eZJQt/roDoKsZnH841ic/3KbXe5OefqdQ4Fi380SQmPOR32drVS8z0OReb996aWqj8leljTOLz9/rmv/5srtmw5u9St5W8fbv56bvPjy/gEPXbDwtO6PYWLfnRGGPMcy/faua8llAhbU6f+aT5/odTzIQpt5W7jWmzLzUvzexbIfEcj5TtS07YtvBcXeg3rwZ02WJlqMjLFuPj47n4Es9NGTv/6ccDA8t/IjNx9HNMv+IaAFL04RVKVart2/+hVavTgx1GtXKsyxZDog+94EgEx5PMAV8yV0pVPk3mFSskErrd+Btu5vj02OX/TkyllKqqQiKhu23Fb0c+Xuf981XplZRSqgoJiYReu3Z6hbc5InFChbeplFKVqdon9A/enIHL5bn6suC1uOXxxMT3KyIkpZQKimqf0HdtWcWp7fyPDldW73bTEzRKqeqr2id0uwUiIz1DuDqO1KqQNh/9/tMKaUcppU6kap/QbeTfhn9F9wfL3U7smPw7FZ8dF388ISmlVFBU+4RuMfkJ/bxLrix3O/O7hebj5pRSJ49qn9ArajialJohMayNUuokVu0TutvtGUjH32h7gXpndv5zM3tv9P+YKqWUquoCSugicrWIrBeRTSIy3M/yS0TkbxFxisjtFR9mydx4EnrBkdzKakVWPd/0VWbzMWoqpVTVVWpCFxErMAO4BugA9BGRohd87wAeAN6t6ABL4zSerpI9u+uUu4253fN3595+A487JqWUCoZAOo67A5uMMVsAROR94CbA9zh3Y8w27zL/jyivRE7vbf8xTcv21G6llAo1gXS5NAN2FphP9paVmYj0FZEkEUlKS6uYR7tJ2BEAduxsUK71pz03xjf96PcfVkhMSikVDIEkdH9DGZZrEHVjzBxjTDdjTLfo6OjyNFGMNdzTd74n7dRyrb82sotv+tlx4yokJqWUCoZAEnoy0KLAfHNgd+WEU3a28EMAmCMRpdQsbvzIOP7XsWVFh6SUUkERSEJfArQTkTYiEgbcBfh/qOMJ9trUcdi9R+jklf1BzFN73eKb7vvrnxUVllJKBUWpCd0Y4wQGAwuAdcA8Y8waEUkQkRsBRORcEUkG7gBmi8iaygz6qOy0HYSFHcbtFm6+qW+Z1u07/fNC8wlxlfMUbqWUOlECuj3SGPMV8FWRstgC00vwdMWcUDaLYA/LxXHYzpnnnBPwerNffo75nfMfNRdz0FkZ4Sml1AlVre8UtRohLOwweUfsZVovrnPh54Y+sLxK9CAppdRxqd4JXWyEhR3myJGwgOonxk6kyU/LC5XdunIbj41KqIzwlFLqhKrWCR1jJ8yey+G80nuOxsRPZfrlvQuVXbItnelDrq+s6JRS6oSq1kMMutx27GG5HHSUfFPRvLdmszi1Lv+99JLiy/7dszLDU0qpE6paJ3SnzY2IIddZE4DXp75Aamo2YmBK79s8lVr2AD+Xmj/23cdweecTGK1SSlWuap3Qjfe2/yO5tUgcNY5Zl1+L6yx/N7YWNuSH+YwYP7ayw1NKqROqWid0S3guAFO7jA6o/hnph7npn694PFGTuVIq9FTrhG6teZgpPFlqvVtW76DlnkWMGD8eOL/yA1NKqSCotgn9gzdnEB52iL/kAr/Lz0/OonVKKo/0Oo0OQ24EbjyxASql1AlWbRP67q2rqH9a4acUXbN+D1ewnq49LqPD5ZcFJzCllAqSapvQaxLO3Kg+vvmUyzt7r1q5puSVlFIqhFXjG4usZNsjARj488/BDUUppaqAapvQncbODmkDQOyYx4McjVJKBV+1TegHa57wx5cqpVSVVi0T+pq//+a39h0BuP2f5aXUVkqpk0O1TOjffzqHJQ3OBiB687IgR6OUUlVDtUzoVpN/cU7ci1OCGIlSSlUdASV0EblaRNaLyCYRGe5nebiIfOBdvlhEWld0oAVtaHUKADfs+LsyN6OUUtVKqQldRKzADDwXeHcA+ohIhyLVHgIyjTGnAi8DEyo60KNih/6LH1ufBUCzdUsrazNKKVXtBHKE3h3YZIzZYozJA94HbipS5ybgLe/0R0BPESl92MNySOvcg33WaB7a/gHxz8+ojE0opVS1FEhCbwbsLDCf7C3zW8cY4wT2Aw2LNiQifUUkSUSS0tLSyhVwq4yNXJezgJpr88q1vlJKhapAbv33d6RtylEHY8wcYA5At27dii0PxPCh3pOgN5RnbaWUCl2BHKEnAy0KzDcHdpdUR0RsQF0goyICVEopFZhAEvoSoJ2ItBGRMOAuYH6ROvOB+73TtwM/GmPKdQSulFKqfErtcjHGOEVkMLAAsAJvGGPWiEgCkGSMmQ+8DswVkU14jszvqsyglVJKFRfQ8LnGmK+Ar4qUxRaYzgXuqNjQlFJKlUW1vFNUKaVUcZrQlVIqRGhCV0qpEKEJXSmlQoQE6+pCEUkDtpdz9ShgXwWGUx3oPp8cdJ9PDsezz62MMdH+FgQtoR8PEUkyxnQLdhwnku7zyUH3+eRQWfusXS5KKRUiNKErpVSIFZJFaQAABIhJREFUqK4JfU6wAwgC3eeTg+7zyaFS9rla9qErpZQqrroeoSullCpCE7pSSoWIapfQS3tgdXUhIi1E5CcRWScia0TkMW95AxH5TkQ2ev+v7y0XEZnq3e+VInJOgbbu99bfKCL3l7TNqkJErCKyTES+8M638T5cfKP3YeNh3vISHz4uIiO85etF5Krg7ElgRKSeiHwkIv943+/zQ/19FpGh3r/r1SLynohEhNr7LCJviMheEVldoKzC3lcR6Soiq7zrTBUJ4LGexphq8w/P8L2bgVOAMGAF0CHYcZVzX2KAc7zTtYENeB7C/QIw3Fs+HJjgnb4W+BrP06HOAxZ7yxsAW7z/1/dO1w/2/pWy708A7wJfeOfnAXd5p18BBninBwKveKfvAj7wTnfwvvfhQBvv34Q12Pt1jP19C3jYOx0G1Avl9xnPIym3AjUKvL8PhNr7DFwCnAOsLlBWYe8r8Bdwvnedr4FrSo0p2C9KGV/A84EFBeZHACOCHVcF7dtnQG9gPRDjLYsB1nunZwN9CtRf713eB5hdoLxQvar2D88Tr36A/2/v/F2jCoI4/hmIKEZEIyhRi3ggtkZSxB+FqEQMoq0gRNR/wEqQVPYiKRQbxUJEQQ0iaSzUOmBAVDTiBUVPoomFP7AKOBY7R16e9y5nPLjbZT6w3Huzm+PNfi/De7N7N+wDxuzD+hXoyGtM+A3+nXbcYeMkr3t2XLs1YLUFN8nZk9WZ+RrDXabbGHAwRZ2BnlxAb4qu1jeZsS8YV9RiS7k0UrA6OuwRsxcYBzao6jSAva63YUW+xzYnI8BZ4LedrwO+aSguDguvv6j4eEw+l4BZ4Lqlma6KSCcJ66yqn4ALwAdgmqDbBGnrXKVZum6y47y9LrEF9IaKUceEiKwC7gFnVPVHvaE1bFrH3naIyGFgRlUnsuYaQ3WRvmh8Jtxx7gCuqGov8IvwKF5E9D5b3vgoIU2yEegEDtUYmpLOi/GvPi7J99gCeiMFq6NBRJYRgvlNVR018xcR6bb+bmDG7EW+xzQnu4EjIvIeuE1Iu4wAayQUF4eF119UfDwmnytARVXH7fwuIcCnrPMB4J2qzqrqHDAK7CJtnas0S9eKHeftdYktoDdSsDoKbMX6GvBaVS9murIFt08QcutV+5CtlvcD3+2R7iEwICJr7c5owGxth6qeU9XNqtpD0O6xqh4HnhCKi8PfPtcqPv4AOGa7I7YAWwkLSG2Hqn4GPorINjPtB16RsM6EVEu/iKy0z3nV52R1ztAUXa3vp4j02xwOZd6rmFYvKixhEWKQsCNkChhu9fX8hx97CI9Qz4Fn1gYJucNHwFt77bLxAlw2v18AfZn3OgWUrZ1stW8N+r+X+V0uJcI/ahm4Ayw3+wo7L1t/KfP3wzYXb2hg9b/Fvm4HnprW9wm7GZLWGTgPTAIvgRuEnSpJ6QzcIqwRzBHuqE83U1egz+ZvCrhEbmG9VvOv/juO4yRCbCkXx3EcpwAP6I7jOIngAd1xHCcRPKA7juMkggd0x3GcRPCA7jiOkwge0B3HcRLhD3MbHVjVEl6bAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "continue_train = True  # set to True to continue to train\n",
    "num_big_steps = 50     # number of small steps\n",
    "num_small_steps = 200  # number of big steps\n",
    "if continue_train:\n",
    "    for k in range(num_big_steps):\n",
    "        for j in range(num_small_steps):\n",
    "            dXY_list = np.append(dXY_list, np.zeros((rep, 1)), axis=1)\n",
    "            dX_list = np.append(dX_list, np.zeros((rep, 1)), axis=1)\n",
    "            dY_list = np.append(dY_list, np.zeros((rep, 1)), axis=1)\n",
    "            for i in range(rep):\n",
    "                minee_list[i].step()\n",
    "                dXY_list[i, -1], dX_list[i, -1], dY_list[i, -1] = minee_list[i].forward()\n",
    "        # To show intermediate works\n",
    "        for i in range(rep):\n",
    "            plt.plot(dXY_list[i, :],label='dXY')\n",
    "            plt.plot(dX_list[i, :],label='dX')\n",
    "            plt.plot(dY_list[i, :],label='dY')\n",
    "            plt.title('Plots of divergence estimates')\n",
    "        display.clear_output(wait=True)\n",
    "        display.display(plt.gcf())\n",
    "    display.clear_output()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save current results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Current results saved.\n"
     ]
    }
   ],
   "source": [
    "overwrite = False  # set to True to overwrite previously stored results\n",
    "if overwrite or not os.path.exists(chkpt_name):\n",
    "    minee_state_list = [minee_list[i].state_dict() for i in range(rep)]\n",
    "    torch.save({\n",
    "        'dXY_list': dXY_list,\n",
    "        'dX_list': dX_list,\n",
    "        'dY_list': dY_list,\n",
    "        'minee_state_list': minee_state_list\n",
    "    }, chkpt_name)\n",
    "    print('Current results saved.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the ground truth mutual information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MI Ground truth is 0.26809677547661703 nats.\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "from scipy.integrate import quad, dblquad\n",
    "from scipy.special import xlogy\n",
    "\n",
    "def fy(y, mu, delta, mix):\n",
    "    return mix*(1/np.sqrt(2*math.pi*(delta**2))*np.exp(-(y-mu)**2/(2*delta**2)))+\\\n",
    "            (1-mix)*(1/np.sqrt(2*math.pi*(delta**2))*np.exp(-(y+mu)**2/(2*delta**2)))\n",
    "\n",
    "lim = np.inf\n",
    "Hy = quad(lambda y: -xlogy(fy(y,mu,rho,p), fy(y,mu,rho,p)), -lim, lim)\n",
    "Hyx = np.log(rho*np.sqrt(2*math.e*math.pi))\n",
    "mi = Hy[0]-Hyx\n",
    "print('MI Ground truth is {0} nats.'.format(mi)) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Apply moving average to smooth out the mutual information estimate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "mi_ma_rate = 0.01            # rate of moving average\n",
    "mi_list = (dXY_list-dX_list-dY_list).copy()    # see also the estimate() member function of MINE\n",
    "for i in range(1,dXY_list.shape[1]):\n",
    "    mi_list[:,i] = (1-mi_ma_rate) * mi_list[:,i-1] + mi_ma_rate * mi_list[:,i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the mutual information estimate after different number of iterations. The red dashed line shows the ground truth, and the green dotted line is the number of iterations where 90% of the ground truth is reached."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeXwV5fX48c/JvpOQsG8JgiigoIAIasUqihv8qlSxbtRWWpcqrbZVa61V29pWv1UrbcWqWDe0aAUs7gpVcQE0IKuEPSYBEpKQfT2/P2YSLyHLTbiTm+W8X6/7ytxZnjlzZzJnZp6ZZ0RVMcYY032FBDsAY4wxwWWJwBhjujlLBMYY081ZIjDGmG7OEoExxnRzlgiMMaabs0TQCBFZLiI/bKd5XScie0WkWESS22OeR8KNc2iw4+gIRORyEXkr2HEEkojMFpEPgzj/Jv8fRGSw2z80WPE1R0Q2iMiUYMfRFt02EYjIThEpczesvSLylIjEtbKMVBFREQlrYwzhwP8BZ6tqnKrmNVH+5w36p4hIpYjsbLA8Z7UljmbiOywhunFuD+R83HkFdQfUFqr6nKqefaTluOt4WCBi6sxa+n9Q1d1u/xp3/HY7YGsk1gUicl+D+Eap6vJgxHOkum0icF2oqnHAicAE4M52nn8fIArY0MJ4sSIy2uf794AdnkVlTAC04QDJ3/+HgOioZxZBoard8gPsBM7y+f5n4DW3eznwQ7c7BCdB7AL2Af8CerjDdgMKFLufSY3MJxJ4CMhyPw+5/Y4GSnymf6+RaVPd4XcCf/bpvxr4FbCzqeVpUE4IcBuwDcgDXgJ6usOigGfd/gXAKpx/yN8BNUC5G9+j7vgKDHO7FwB/A153x/kI6OsuYz6wGTjBJ466GIqAjcB33P7HuvOpccsp8PntHnB/573AP4Bod1gK8Job8wHgAyCkieV/GNgDHATWAKf5DIsGnnbj3QT8AshsKWZ32GzgQ5/vCvwY2OqWNw8Qd9gwYAVQCOQCL7r9/+dOV+Iu+6WNxD8b+ND9LfJxDgLObWZbvht4tsE29H33N8h3Y5wArHN/v0cbzOsj4K9urJuBM32G9wCeALKBr4H7gNAG0/7FXSf3efT/EEbT2+cxwNvu/LcAl/hMvwD4O7DMnddZwPnAF+62sQe4u8E8TwVWur/THncZ5wBVQKU776UN10NTy+kOmwJkArfg7FOyge/7zPM8nG2tyP2Nb/V8f9heO96O9mmw0gbhHIXc635fzjeJ4BogAxgKxAGvAM803DCbmc89wCdAb6CXu1Hd68/0PsNT3Y0wFGenucXdiP1NBHPdGAa6G+hjwAvusB8BS4EYt/xxQELD38GnrIaJINedJgp4D2cndZVb1n3A+z7Tfhfoj5OYLsX5Z+znDpuNz07V7fcQsAToCcS7cf7BHfYHnMQQ7n5Ow93pNrL8VwDJODuQW4AcIModdj/ODjrJ/X3WcWgi8Dtm97d5DUgEBgP7gWnusBdwkneI+1ud2thv2kT8s3F2PNe6v+t1ODuXuiRzyLqn8UTwD3e+Z+PsPF/F2SYH4OyMTveZVzXwU/d3vRQnIdQdOLyKs/3EutN/BvyowbQ/cX/raI/+H8Ia2z7dmPbgJL0wnDP9XGCUz/ZaCJzisx6mAMe534/HOeD4f+74g3F2xpe5v0UyMNanrPsaxFe/HlpYzinu73SPW+55QCmQ5A7Pxj1YwdkuT/R8f+j1DDrqx11pxTiZfhfOkW3d0Wb9Bga8C1zvM90InH/KsJY2XHf8bcB5Pt/Pwd2Bt2bDB95xp70fZ4fSmkSwiUOP6vr5LMM17kZ6fCPTHfKP5vZrmAge9xn2E2CTz/fjcI/um4grHZjhds/m0J2q4Ox0j/LpNwnY4XbfAyymmR1oM/PNB8a43duBc3yG/RCfRNDKmJVDd/AvAbe53f8C5gMDGynTn0SQ4fM9xp2mb2PrnsYTwQCf4Xn4nHkALwNzfeZVn2Tcfp8BV+KcKVbgs4PH2Um+7zPt7hZ++4D8PzS2feIkrQ8aTPMY8Buf7fVfLcT3EPAXt/t24D9NjLeA5hNBc8s5BSjzXU6cZHyy270b5wAtobXbdls/3b2O4P+paqKqDlHV61W1rJFx+uMkijq7cHagffycR2PT929DrP/C+Ue7DOdSTmsMAf4jIgUiUoCTGGpwluEZ4E1goYhkicif3Eo7f+316S5r5Ht9BbyIXCUi6T5xjMa5xNOYXjg7vDU+47/h9gfnUl4G8JaIbBeR25oKUERuEZFNIlLoltPDZ779cY4i6+xpMG1rYgbnbKNOqc/y/wInuX3m3l1yTTNlNFuuqpa6na25ucHv9QR8re4eyVW3zQ7BOYLN9vk9HsM56q1zyO/XiED9PzRmCDCxLjY3vstxLlc2Gp+ITBSR90Vkv4gU4lw2q1u/g3B26G3R0nLmqWq1z3ffbeVinLOEXSKyQkQmtTEGv3X3ROCPLJwNrM5gnNO6vThHJ22ZPqsNcbyMcz1zu6ruamnkBvbgXFNO9PlEqerXqlqlqr9V1ZHAZOACnEs74N/y+UVEhgCPAzcCyaqaCKzH2Tk2Nq9cnB3UKJ+Ye6hTuY+qFqnqLao6FLgQ+JmInNnIfE8DfglcgnPqnYhzeaBuvtk4l4TqDGpFzH5T1RxVvVZV++Mc7f0tgHcKleAkzTp9mxrRTwNExHcZ67bZPThnBCk+6yRBVUf5jNvSNhOo/4fG5rUHWNFgO49T1euameZ5nMuPg1S1B84lNPEp7yg/591Qm5dTVVep6gycBPsqzpmlpywRtOwF4KcikubeXvp7nIq+apxrwLU49QfNTX+niPQSkRTgLlp/RI+qlgDfxrl00Vr/AH7n7thwY5nhdp8hIse5d1AcxLlkVONOt5fml601YnH+efa78/0+ztF1nb3AQBGJAFDVWpyd8F9EpLc7zQAROcftvkBEhrk7rINuzDUcLh4nce8HwkTkLiDBZ/hLwO0ikiQiA3B2+v7G7DcR+a6I1CWcfLfcQP3O6cAsEQkXkfHAzCMoC5wd0E1ued/FqZdapqrZwFvAgyKSICIhInKUiJzeirID8v/gavi7vQYcLSJXurGHi8gEETm2mTLigQOqWi4iJ+HckVfnOeAsEblERMJEJFlExjYx74batJwiEiHO8yk9VLWKb7ZtT1kiaNmTOJdP/odTEVqOcy287hT9d8BH7qnoyY1Mfx/OXT7rgC+Bz91+raaqq1W1LaeqD+Mc9bwlIkU4lVgT3WF9gUU4G9wmnIrTZ32mmyki+SLySFti9ol9I/Ag8DHOP9FxOHeY1HkPp8I+R0Ry3X6/xLn884mIHMSpJxnhDhvufi92y/ybNn4P95s4dzV9hXN6Xs6hlwfuwbmDY4db3iKco15/Ym6NCcCnIlKMsy5uVtUd7rC7gafdbeiSNpT9a5wj13zgtzhHuUfiU5zfNxdn+56p39zTfxUQgXNXSz7O79WvFWUH7P+BBtunqhbhVIbPwjn6zgH+iHODRFOuB+5x/y/uwufoW1V341yiuQXnLqR0YIw7+AlgpLvOXg3wcl4J7HS3+R/j3Ozgqbq7DowxOE+2ArNUtTVHucZ0anZGYLo1EeknIqe4lzlG4Bz9/SfYcRnTntrUNIIxXUgEzp0vaTi3Ei/EuZXYmG7DLg0ZY0w3Z5eGjDGmm+t0l4ZSUlI0NTU12GEYY0ynsmbNmlxV7dXYsE6XCFJTU1m9enWww2jWnkLn7sRBPQa1MKYxxrQPEWnyQdROlwg6gyv/cyUAy2cvD24gxhjjB0sEHrjzW+39WgNjjGk7SwQeOGtoQF8UZowxnrK7hjywPX872/MD/jZHY4zxhJ0ReOCaxU4Lw1ZHYIzpDCwReOC3U34b7BCMMcZvlgg8cHqqtVdmjOk8LBF4YEvuFgBGpIxoYUxjTE2tIgIHSyvJKyqnvKqGvQVl9EuK4WBZJZXVtUw6ug9hoYGt0iwsrSQsVIiNPPSFfKqK73t5VJWqmlpCQ4TC0koSoiMICREqq2sJESitqKZHTMRh0xSWVlKrSlS4s5s9WFpJQWkFeUUV1NYqhWWVjE1NJq+ogn2FZUSEhdCrRzRf55WgKCMHJpESH0VURBg1tbWUlFcTERZCVETgd9uWCDzwo9d+BFgdQWdQWV1DaIiQX1xJRHgICdER9f/0IT7/2CUVVRSXVVFTq0RFhKIK8dHhhIgQGiKUVlaTX1xBYWklqlBR5bxLpKKqhrQ+CVRV11BWVUNyXBRhoUJNrdIjJoKisioKSyvZklVASUU1AkSEhZBfXMGu3GJOSEshLFQoLq+mrKIaEUjtHU9UeBgicKConPySSo7u34NByXFs23uQ6IhQisqqWJWxj75JMVRV11JWWcPB0krKKqspKK2kplYJC5H6HVhcVBgJ0RGEhYawOauA4wf3pKKqhsKyShJjIlFVPv5qLxk5Bzl2YCI9YiIJERjRP5EBPWNJiIlgaJ8Eduw7yN6CMkorq8nOLyVEoKS8mhpVDpZWEhcVTl5xOcVlVZRXOuVXVtVQUV3b4rqKiQwjLiqcXglRxEaFc6ConKiIMJJiI4kIC2FPbjF5xeUMSokjI9t5CV2vhCgAqqprKamoJi4qjIrqWmIjw9i1v4hahaTYSEorqg6JoUdMBCEi9O4RTVZ+CUVlVY3GFCJQqxAVHspRfRMor6yhqqaW3bnFbd8oG0iIDqeorAoFbj7/OM47cXDAyq7T6RqdGz9+vHb0J4tX7lkJwORBk4McScemquwtKGPT1/kcLK2kpKKa0JAQhvaJZ2xaCgUlFcRHhRMVEcbBskpiI8Mor6ph9/5iPtu6j41f55McF0VBibMD7pcUS0pClLMDK60kKjyUiqoaCkorEaCqxvlHr6iqoaZWOVhWSWFJZZPvHOzTI5rK6lrCw0LYV9jY66w7DwFio8IBpW9iDFERYazffaBV09f9TnFRYUSFh1FUVtnsDjw8NIS6XCpQP26vhCjS+iQQGRYKQExkKOGhIQxKiSM6Iozw0BDySyqIiwqnb2IMu3OLKSipoKisij15xVRV15JbVE5YiKBAZl4J4Oy8eyVEUV2jlFRUsf9gOX0TowFQhZ5xkURFhBESItTWKqm94zlQVI4CYSFCaUU1g1LiWLszj5yCUgb0jCUmKpz4qHAGJcdSXlVDj5gIyiprqK6p5YPN2Rw/JJnYyDD2FpRRWFpJcXkVsVHh7Mkt5ruThyII5VU1RIaFUFZZw+BecSTFRlJdU0tWfglVNd8k5MTYSEoqqugRE0FEWCgZOYXs3FfE3oIyBiY72/bY1BRSe8f7vd4OWYcia1R1fKPDLBGYQKnbliqqa/l4Sw4FJZW8mb6HqPBQ98ivloSYCEJDhOqaWgpLKzlQXNFiuZHuDr2h8NAQIsNDiAoPY0ByLF/uyqNWISYijMrqGhJiIjhQXEHPuEj6uDsEQep3ToNT4kiOj6K8qobeCVEcKK5g+74i+iXGUFNbS2ZeCXlF5STHRxERFsL4Yb0JCxGy80upqK6hoqqGEBHCw0JIio0kJSGKyLBQFGenk1dUTkbOQRJjI4iJCCM6Iow8nx1PVn5p/T/9iUNTSHFjKSqrJCkukviocA4UVzjvtKxVQsQ5Ks4pKCO/uIKYyDAiwkLoGRfFrv1F7M4tZlByHCLUn7EM6BlLSIgQGxnW5KWVWvdovbC0ktKKahJjnaP9ODcJ1501Rbg7bt/1XVJRzY69B8krruCL7bmcMDSFmIgwkuIiGdongdAQOWxeIXJov0BQVWqVw+ZnvhG0RCAi03BeJxcK/FNV728w/C/AGe7XGKC3+4LwJnWGRLB+33oARvdu0+ttO43K6hq2fF1AVn4pn27dx0ebcwgPDaFWlZrab7ar+OhwRg1MolYVRCgsca/NRoVzVJ8ETj66NwnREZRWVJMUF8lXWYVsySpgX2EZ4WEhVFXX0i8phvKqGiLCQujdI5rxR/UiJT7qkOuyjWl4vdeY7qq5ROBZHYH7MvR5wFScd8KuEpEl7ntgAVDVn/qM/xPgBK/iaU83LnPef97V6ggqqmrYklXAlqwCPtqcw/acg/Wn+6EhwvB+PTh2YCJR4WGMG5rCoJQ4esREtLqSb9KIKCaN6BOQmC0JGNMyLyuLTwIyVHU7gIgsBGbgvPS6MZcBv/Ewnnbz56l/DnYIAVNeWc2Kjdms3JzDmu259dfZQwTOHjOIcUf1onePaIb2iT/s0oExpnPwMhEMAPb4fM8EJjY2oogMwXlV4HtNDJ8DzAEYPDjwNeaBNmHAhGCHcMQqq2t4Mz2TZ1Z8RWFpJT3jIjl/3GDGDElm5KAkEmMjgx2iMSZAvEwEjZ2TN1UhMQtYpKqH1wgCqjofmA9OHUFgwvNOek46AGP7jg1yJG3z+fZc/rw4nQPFFYwalMRN541m0oi+VhFnTBflZSLIBHzfzDIQyGpi3FnADR7G0q7mvjEX6Jx1BK9/sZtHl62nf89YfvH/xjI2NdmusxvTxXmZCFYBw0UkDfgaZ2f/vYYjicgIIAn42MNY2tVD0x4KdgitVlOrPPr6epZ9vptxQ1O44+ITiYsKb3lCY0yn51kiUNVqEbkReBPn9tEnVXWDiNwDrFbVJe6olwELtbM90NCMznZJSPWbJHDxyWn84Mxj7TKQMd2Ip01MqOoyYFmDfnc1+H63lzEEw6qvVwGdp9L4xY+2sezz3Vx0chpzpo4MdjjGmHZmbQ154Odv/xzoHHUEiz7ezlPvb+GUY/py7VnHBjscY0wQWCLwwKPnPRrsEPyyeNVOHn9nEycN782vLj7Rk0f/jTEdnyUCD3SGpiW27z3IY29tZKKbBKxOwJjuy95Z7IGVe1bWt0DaEVVW1/CnV9OJjw7n1uljiAy3J4KN6c7sjMADd7x7B9Bx6wie+99Wduwr4p5Z40mIiQh2OMaYILNE4IHHLngs2CE0KaeglJc/2cG3R/dn4vDANOxmjOncLBF4oCO/ovKJdzcTInDNmccEOxRjTAdhdQQeWLFzBSt2rgh2GIfZlJnP/zZm893JR9ErITrY4RhjOgg7I/DAb5Y7rWl3tDqCFz/aRnx0ODMnDQ12KMaYDsQSgQeenPFksEM4zJ7cYj75ai+XnTaM6Ahb7caYb9gewQNDkzreEfcrn+4gLDSE6eNTgx2KMaaDsToCD7yz/R3e2f5OsMOol19cwdtrM5k6ZiBJcfZCGWPMoeyMwAP3/e8+AM4aelaQI3EsWbWT6ppaLj45LdihGGM6IEsEHnjmO88EO4R65ZXVLF2zi0kj+jAwOS7Y4RhjOiBLBB4Y1GNQyyO1k/fWZ1FUVsXFJ3e8egtjTMdgdQQeeCPjDd7IeCPYYaCqLFm1k6F9Ehg1KCnY4RhjOig7I/DA/R/eD8C0YdOCGsf6Pfns2FfE3AuOs/cOG2OaZInAAwtnLgx2CAC8tnoXsZFhnDGqf7BDMcZ0YJ5eGhKRaSKyRUQyROS2Jsa5REQ2isgGEXney3jaS9+4vvSN6xvUGIrLq/hocw5nHT+QKHuAzBjTDM/2ECISCswDpgKZwCoRWaKqG33GGQ7cDpyiqvki0tureNrT0i1LAbhwxIVBi+GjzTlU1dRy5vEDghaDMaZz8PJQ8SQgQ1W3A4jIQmAGsNFnnGuBeaqaD6Cq+zyMp908+PGDQHATwYqN2fRLiuHofj2CFoMxpnPwMhEMAPb4fM8EJjYY52gAEfkICAXuVtXDbrcRkTnAHIDBgwd7EmwgLbpkUVDnX1xeRfqOXC6amGaVxMaYFnmZCBrbA2kj8x8OTAEGAh+IyGhVLThkItX5wHyA8ePHNyyjw0mJSQnq/D/fnktNrTJphL14xhjTMi8rizMB3yerBgJZjYyzWFWrVHUHsAUnMXRqr2x6hVc2vRK0+X+WsY+4qHCOGZAYtBiMMZ2Hl4lgFTBcRNJEJAKYBSxpMM6rwBkAIpKCc6lou4cxtYtHPn2ERz59JCjzVlXSd+RyQloyoSH2vKAxpmWeXRpS1WoRuRF4E+f6/5OqukFE7gFWq+oSd9jZIrIRqAF+rqp5XsXUXhbPWhy0eWfnl7L/YDmXnpIctBiMMZ2LpzeYq+oyYFmDfnf5dCvwM/fTZfSICt6dOut2OXl0zBBLBMYY/9i1Aw+8uP5FXlz/YlDmvW7XARJjIxiUYi2NGmP8Y4+ceuDvq/8OwKWjL23X+aoqa3flcfyQZLtt1BjjN0sEHlh2+bKWR/JAdn4puQfLOd4uCxljWsESgQdiwmOCMt9v6gd6BmX+xpjOyeoIPPDsumd5dt2z7T7ftTvzSIqNtPoBY0yr2BmBB/75+T8BuOL4K9ptnqrKul0HOH5IT6sfMMa0iiUCD7x95dvtPs+s/FJyi8o5PtXqB4wxrWOJwAPhoeHtPs+6+gGrKDbGtJbVEXhgQfoCFqQvaNd5bv66gPjocAYlx7brfI0xnZ8lAg8EIxF8lVXI0f16WP2AMabV7NKQB5bPXt6u8yuvqmHnviImnnJUu87XGNM12BlBF7Atp5BaVY7ub28jM8a0niUCDzy+5nEeX/N4u83vq6xCAEb0t/cPGGNazxKBB17c8CIvbmi/Rue+yiogOT6S5PiodpunMabrsDoCD7xz1TvtOj+notjOBowxbWNnBJ1cSXkVmQdKrH7AGNNmlgg88LdVf+Nvq/7WLvPamu3UDxxt9QPGmDayROCBpV8tZelXS9tlXlvciuKj+9kZgTGmbTxNBCIyTUS2iEiGiNzWyPDZIrJfRNLdzw+9jKe9vH7567x++evtMq+vsgrolxRDQkxEu8zPGNP1eFZZLCKhwDxgKpAJrBKRJaq6scGoL6rqjV7F0dV9lV3IsQPsspAxpu28PCM4CchQ1e2qWgksBGZ4OL8O4+FPHubhTx72fD4FJRXsKyyz+gFjzBHxMhEMAPb4fM90+zV0sYisE5FFIjLIw3jazbs73uXdHe96Pp8tWQUAjLA7howxR8DL5wgaa/1MG3xfCrygqhUi8mPgaeDbhxUkMgeYAzB48OBAxxlwSy5b0i7z+SqrEAGO6muJwBjTdl6eEWQCvkf4A4Es3xFUNU9VK9yvjwPjGitIVeer6nhVHd+rVy9Pgu2MvsoqYFBKHDGR9lygMabtvEwEq4DhIpImIhHALOCQQ2UR6efzdTqwycN42s0DKx/ggZUPeDoPVWVLVqG1L2SMOWKeHUqqarWI3Ai8CYQCT6rqBhG5B1itqkuAm0RkOlANHABmexVPe/o482PP57GvsIzC0kp7otgYc8Q8vaagqsuAZQ363eXTfTtwu5cxBMPLl7zs+Ty+sieKjTEBYk8Wd1JfZRUSFiIM7RMf7FCMMZ2cJQIP3P/h/dz/4f2ezmNrdiGpveOJCAv1dD7GmK7PbjfxQHpOuqflqyoZOYWcckxfT+djjOkeLBF4YOHMhZ6Wv6+wjKKyKobZ8wPGmACwS0OdUEbOQQCG9U0IciTGmK7AEoEH7l1xL/euuNez8jOyCwkRYWgfSwTGmCPXYiIQkaNF5F0RWe9+P15E7vQ+tM5rS94WtuRt8az8jJxCBqfEERluFcXGmCPnTx3B48DPgccAVHWdiDwP3OdlYJ3Zsxc962n5GTkHOXFoiqfzMMZ0H/5cGopR1c8a9Kv2IhjTsryicg4UV1hFsTEmYPxJBLkichRuy6EiMhPI9jSqTu6u9+/irvfvannENsjIcZ4oHmavpjTGBIg/l4ZuAOYDx4jI18AO4HJPo+rk9hzc0/JIbZSRfdBpetoqio0xAeJPIlBVPUtEYoEQVS0SkTSvA+vMnprxlGdlZ+QUMiA51pqeNsYEjD+Xhl4GUNUSVS1y+y3yLiTTnIycg1Y/YIwJqCYPK0XkGGAU0ENELvIZlABEeR1YZ3b7O06Dqn846w8BLbewtJJ9hWVMnzAkoOUaY7q35q4vjAAuABKBC336FwHXehlUZ5dXludJuXUVxcPtjMAYE0BNJgJVXQwsFpFJqur9m1a6kPkXzvek3Ixsp2kJe0exMSaQ/Klx/EJEbsC5TFR/SUhVr/EsKtOojJxC+iZGEx8dHuxQjDFdiD+Vxc8AfYFzgBU4L6EvanaKbu7Wt27l1rduDXi5GTmFVlFsjAk4fxLBMFX9NVCiqk8D5wPHeRtW51ZWVUZZVVlAyywpryLrQKk9SGaMCTh/Lg1VuX8LRGQ0kAOk+lO4iEwDHsZ5ef0/VbXR13a5Tyv/G5igqqv9Kbsjm3f+vICXaU1PG2O84s8ZwXwRSQJ+DSwBNgJ/amkiEQkF5gHnAiOBy0RkZCPjxQM3AZ+2Iu5up/6OITsjMMYEWItnBKr6T7dzBTC0FWWfBGSo6nYAEVkIzMBJJL7uxUksgb+oHiRz35gLwEPTHgpYmRnZhaQkRJEYGxmwMo0xBvxIBCKSCFyFczmofnxVvamFSQcAvo3uZAITG5R9AjBIVV8TkSYTgYjMAeYADB48uKWQuyR7otgY4xV/6giWAZ8AXwK1rShbGumn9QNFQoC/ALNbKkhV5+M0fMf48eO1hdGDLpBnAgClFdXsyS3m9JH9AlquMcaAf4kgSlV/1oayM4FBPt8HAlk+3+OB0cByEQHnFtUlIjK9K1QYB1JGTiEKjBiQGOxQjDFdkF/PEYjItSLST0R61n38mG4VMFxE0kQkApiFU9kMgKoWqmqKqqaqairOWUeXSAI3/PcGbvjvDQErb8vXBQAc3d8SgTEm8Pw5I6gE/gz8im8u7SgtVByrarWI3Ai8iXP76JOqukFE7gFWq+qS5qbvzKLDowNa3uavC+ibGE2PmIiAlmuMMeBfIvgZzkNlua0tXFWX4dQx+PZr9NVdqjqlteV3VA+c/UDAylJVNmbmMyY1OWBlGmOML38uDW0ASr0OxDRuX2EZB4orGDkwKdihGGO6KH/OCGqAdBF5H6io6+nH7aPd1pylc4DAtEK6KdOpHzjWEoExxiP+JIJX3Y/xU3J04C7jbMzMJzI8lKF94gNWpjHG+PLnyeKn2yOQriSQbybblJnPiP49CA3x5yqeMca0XnOvqnxJVS8RkS/xeRCsjqoe72lkhpLyKjJyCpl1yrBgh2KM6cKaOyO42f17QXsE0pV8f/H3AXhqxlNHVM6Xuw9Qq+mXtkAAACAASURBVDAmze4YMsZ4p8nrDaqa7XZer6q7fD/A9e0TXuc0KGEQgxIGtTxiC9buzCM8NMTuGDLGeMqfyuKpwC8b9Du3kX7Gdc8Z9wSknPSdeYwclEREWGhAyjPGmMY0eUYgIte59QPHiMg6n88OYF37hdg9HSytZPveg4y1B8mMMR5r7ozgeeB14A/AbT79i1T1gKdRdXJXvHIFAM9e9Gyby1i7Kw/Anig2xniuyUSgqoVAoYjcCeSoaoWITAGOF5F/qWpBewXZ2YxIHnHEZazdmUdUeCgjrKE5Y4zH/KkjeBkYLyLDgCdwWhB9HjjPy8A6s1+f/usjLmPtzjxGD+5JWKg9P2CM8ZY/e5laVa0GLgIeUtWfAvaGFA/lFZWzO7fY6geMMe3Cn0RQJSKX4byu8jW3X7h3IXV+sxbNYtaiWW2efuWWvQBMGNY7UCEZY0yT/Lk09H3gx8DvVHWHiKQBba8F7QbG9h17RNN/uDmbgT1jGdIrLkARGWNM0/xpa2ijiPwSGOx+3wHc73Vgndltp97W8khNOFhaybqdB/ju5KG4r/A0xhhPtXhpSEQuBNKBN9zvY0Wky75dLNg+/movtaqcdqxVwxhj2oc/dQR3AycBBQCqmg6keRhTp3fxSxdz8UsXt2naDzfn0KdHNMP6JgQ4KmOMaZw/iaDafabA12GtkTZGRKaJyBYRyRCRw66XiMiPReRLEUkXkQ9FZKQ/5XZ0kwZOYtLASa2erqCkgjXb9nPqsX3tspAxpt34U1m8XkS+B4SKyHDgJmBlSxOJSCgwD6etokxglYgsUdWNPqM9r6r/cMefDvwfMK2Vy9Dh3Dr51jZN996XX1NTq5w95sgbrDPGGH/5c0bwE2AUzmsqnwcKgbl+THcSkKGq21W1ElgIzPAdQVUP+nyNxc8zja5IVXlrbSZH9+9Bam97G5kxpv34c9dQKfAr99MaA4A9Pt8zgYkNRxKRG4CfARHAt1s5jw5p+gvTAVhymf916l/uPsCOfUXcfP5xXoVljDGN8ufSUFs1dpG7sTedzQPmuZef7gSuPqwgkTnAHIDBgwcHOMzAOzPtzFZPs2TVTuKjw/n2cQM8iMgYY5rmZSLIBHwvdg8EspoZfyHw98YGqOp8YD7A+PHjO/zlo5tPvrnlkXzsKyzjo817uejkNKLC7d0Dxpj25WWLZquA4SKSJiIRwCycBuvquZXPdc4HtnoYT4f10spthAhMHz8k2KEYY7qh5l5e/1eaqbxV1ZuaK1hVq0XkRuBNIBR4UlU3iMg9wGpVXQLcKCJnAVVAPo1cFuqMzn3uXABev/z1FsfNPVjOG1/sYeqYgfRJjPE6NGOMOUxzl4ZWH2nhqroMWNag310+3a27htJJXHj0hX6P+9wHW6lVZdYpwzyMyBhjmtbci2mebs9AupLrJ1zv13gZ2YW88cVupk9IpW+SnQ0YY4KjuUtDzd77qKrTAx9O91FRVcOfFqeTGBvJ5d8a3vIExhjjkeYuDU3CeQ7gBeBTGr8dtP1t2QJTphza75JL4PrrobQUzmvkxWmzZzuf3FyYOfPw4dddB5deCnv2wJVXHj78llvgwgudef/oR4cPv/NOOOssSE+HuXNZu3ctAGP6jHGG//73MHkyrFwJd9wBQO6BEm4oKmdo7wTiz5gHY8fCO+/AffcdXv5jj8GIEbB0KTz44OHDn3kGBg2CF1+Evzdy49WiRZCSAgsWOJ+Gli2DmBj429/gpZcOH758ufP3gQfgtdcOHRYdDa+7dSH33gvvvnvo8ORkePllp/v22+Hjjw8dPnAgPOu2aj53rvMb+jr6aJg/3+meMwe++urQ4WPHwkMPOd1XXAGZmYcOnzQJ/vAHp/viiyEv79DhZ54Jv3bfKHfuuVBWdujwCy6AW90nxRtud9Dhtr3DNLLtHeKhh2zbg+657floLhH0xWke4jLge8B/gRdUdUOzJRp6xTT/Qpn8kgpyi8pJiY8iPtre8WOMCS5Rbfm2fBGJxEkIfwbuUdW/eh1YU8aPH6+rVx9xPXbQrMrYx29eXM3owT2577IJRITZcwPGGO+JyBpVHd/YsGYfKHMTwPk4SSAVeAR4JdABdhcbM/O5999rSOsdz28uGWdJwBjTITRXWfw0MBp4Hfitqq5vt6g6uSkLpgCwfPby+n6fbd3HH1/9gp7xUdx32UnERtolIWNMx9DcGcGVQAlwNHCTT/v4Aqiq2ptTmjB77Oz67rLKap7731b+/fF20nrH89tLx5MUFxm84IwxpoHmniPwsvmJLm322NmoKis2ZPGPtzZyoLiCaScM4vpzRhFpbQkZYzoYLxud67Y2f53Hk+9vZu2OAob2SeDX3x3HyIFJwQ7LGGMaZYkgQKpravlocw6LV+1kwa4fESLCk+cu4fxxQwgN6RiPYBhjTGMsERyhr/NKeP2L3by3/mvyiirolxTDlcd9n9GDezJ9fGqwwzPGmBZZImiDmtpaPv1qH699vpvPt+0nNEQ4YWgKN58/hPFH9SY05Ixgh2iMMX6zRNAKVTW1/OfTHSxdvYt9hWX0SojislOHceGEIfSMi6ofr7SqFICYcGtIzhjT8Vki8FNGdiEPLFnLjn1FjElNZs5ZxzL5mD6Ehhx+c9V5zzntfvg+R2CMMR2VJQI/rNm2n7tfWk1cVDi/vXQ8Jx/dp9nxrxt/XTtFZowxR84SQQvSd+Ry18JVDOkVz+8vP4nE2JYfBrt09KXtEJkxxgSGJYJmFJRU8MdX0+mXFMMfrzzZ75ZCC8sLAegR1cPL8IwxJiAsETTj729upKisit9976RWNRc9Y+EMwOoIjDGdg6eJQESmAQ/jvLz+n6p6f4PhPwN+CFQD+4FrVHWXlzH5a2NmPss3ZPG904YxtE/rmlW6aeJNHkVljDGB51kiEJFQYB7Oy20ygVUiskRVN/qM9gUwXlVLReQ64E9Ah7jA/sIHW0mMjeCSyUe1etqLjr3Ig4iMMcYbXjYsdxKQoarbVbUSWAjM8B1BVd9X1VL36yfAQA/j8dvXeSV8lrGfC8cNITqi9bkytzSX3NJcDyIzxpjA8zIRDMB553GdTLdfU36A8+6Dw4jIHBFZLSKr9+/fH8AQG/dG+h5CQ4Tzxg1u0/QzX5rJzJeaf0eoMcZ0FF7WETTW0lqj78UUkSuA8cDpjQ1X1fnAfHBeVRmoAJuYFys2ZjFuaMohTwu3xi2TbglwVMYY4x0vE0EmMMjn+0Agq+FIInIW8CvgdFWt8DAev+zJLWZvQRmzThnW5jIuHHFhACMyxhhveXlpaBUwXETSRCQCmAUs8R1BRE4AHgOmq+o+D2Px2+c7nGv7Jw5NaXMZOcU55BTnBCokY4zxlGdnBKpaLSI3Am/i3D76pKpuEJF7gNWqugT4MxAH/Nt9FeZuVZ3uVUz++GJ7Lv17xtA3se0Nxs1aNAuw5wiMMZ2Dp88RqOoyYFmDfnf5dJ/l5fxbq7qmlnW7DnDGcf2PqJzbTr0tQBEZY4z37MliH9v2HqS0spqxqW2/LAQwbdi0AEVkjDHesxfU+8jIdtoIOrr/kbURtKdwD3sK97Q8ojHGdAB2RuAjI+cgcVHh9OkRfUTlXPmfKwGrIzDGdA6WCHxk5BQyrG8CbsV1m935rTsDFJExxnjPEoFLVdm1v5hzTxjU8sgtOGtoh6oDN8aYZlkdgaugpJKKqhr6Jx35e4a3529ne/72AERljDHeszMCV06B0/Zd3wAkgmsWXwNYHYExpnOwRODKzncSQb8jeJCszm+n/PaIyzDGmPZiicBVd0bQJwCJ4PTURtvOM8aYDsnqCFzZ+aUkx0cSGR56xGVtyd3CltwtAYjKGGO8Z2cErpyC0iNqX8jXj177EWB1BMaYzsESgSs7v5QxqckBKev3Z/4+IOUYY0x7sEQAVFbXkHuwPCAVxQCTB00OSDnGGNMerI4A2FdYhhKYimKA9fvWs37f+oCUZYwxXrMzAiCnoAyAfgF4hgDgxmU3AlZHYIzpHCwR4PMMQYASwZ+n/jkg5RhjTHuwRIBzx1BEWAhJcZEBKW/CgAkBKccYY9qD1RHgnBH0TYwh5AhbHa2TnpNOek56QMoyxhiveZoIRGSaiGwRkQwROez9jSLyLRH5XESqRWSml7E0Jye/NCBtDNWZ+8Zc5r4xN2DlGWOMlzy7NCQiocA8YCqQCawSkSWqutFntN3AbOBWr+JoiaqSXVDK6ME9A1bmQ9MeClhZxhjjNS/rCE4CMlR1O4CILARmAPWJQFV3usNqPYyjWUVlVZRWVAf0jGBs37EBK8sYY7zmZSIYAPi+uDcTmOjh/NokuyBwrY7WWfX1KsAqjU3XV1VVRWZmJuXl5cEOxbiioqIYOHAg4eHhfk/jZSJorOZV21SQyBxgDsDgwYOPJKbD5Li3jvZNPLL3FPv6+ds/B+w5AtP1ZWZmEh8fT2pq6hG/4tUcOVUlLy+PzMxM0tLS/J7Oy0SQCfi+93EgkNWWglR1PjAfYPz48W1KJk0J5Atp6jx63qMBK8uYjqy8vNySQAciIiQnJ7N///5WTedlIlgFDBeRNOBrYBbwPQ/n1ybZ+aUkxkYQHRG4n2J079EBK8uYjs6SQMfSlvXh2e2jqloN3Ai8CWwCXlLVDSJyj4hMBxCRCSKSCXwXeExENngVT1OyC0oDWj8AsHLPSlbuWRnQMo0xxiuePkegqstU9WhVPUpVf+f2u0tVl7jdq1R1oKrGqmqyqo7yMp7GBPoZAoA73r2DO969I6BlGmMa9/DDDzN69GhGjRrFQw99c+v2gQMHmDp1KsOHD2fq1Knk5+cD8PLLLzNq1ChOO+008vLyANi2bRuzZs1q17h37tzJ6NFHdvUgLi4uILF06yeLq2tq2VcYuOan6zx2wWM8dsFjAS3TGHO49evX8/jjj/PZZ5+xdu1aXnvtNbZu3QrA/fffz5lnnsnWrVs588wzuf/++wF48MEH+eSTT7jqqqt4/vnnAbjzzju59957/ZpnTU2NNwsTRN06Eew/WE6tasDPCEakjGBEyoiAlmlMZzBlwRQWpC8AoKqmiikLpvDsumcBKK0qZcqCKby4/kUACssLmbJgCq9segWA3NJcpiyYwtItSwHIKc5pcX6bNm3i5JNPJiYmhrCwME4//XT+85//ALB48WKuvvpqAK6++mpeffVVAEJCQqioqKC0tJTw8HA++OAD+vXrx/Dhw5ucT1xcHHfddRcTJ07k448/Zs2aNZx++umMGzeOc845h+zsbAAef/xxJkyYwJgxY7j44ospLXVuRtm7dy/f+c53GDNmDGPGjGHlSufScU1NDddeey2jRo3i7LPPpqzMaQl527ZtTJs2jXHjxnHaaaexefNmAHbs2MGkSZOYMGECv/71r/1ZJX7p1okg0K2O1lmxcwUrdq4IaJnGmMONHj2a//3vf+Tl5VFaWsqyZcvYs8d5fGnv3r3069cPgH79+rFv3z4AfvOb33DOOefwzjvvcNlll3Hfffe1uFMtKSlh9OjRfPrpp0ycOJGf/OQnLFq0iDVr1nDNNdfwq1/9CoCLLrqIVatWsXbtWo499lieeOIJAG666SZOP/101q5dy+eff86oUc5V8K1bt3LDDTewYcMGEhMTefnllwGYM2cOf/3rX1mzZg0PPPAA119/PQA333wz1113HatWraJv374B+x1FNaB3Y3pu/Pjxunr16oCUtezz3Tz83y955qZv07tH4J4jmLJgCmDPEZiub9OmTRx77LFBjeGJJ55g3rx5xMXFMXLkSKKjo/nLX/5CYmIiBQUF9eMlJSXV1xPUefrppykoKGDixIk88MADJCUl8fDDDxMTc+jBYVhYGBUVFYSGhrJ+/XomT57M0KFDAeeovl+/frz11lusWLGCO++8k4KCAoqLiznnnHP4xz/+Qa9evcjMzCQy8psWjnfu3MnUqVPrL2X98Y9/pKqqirlz59KrVy9GjPjmqkJFRQWbNm0iOTmZnJwcwsPDOXjwIP3796e4uPiw36Sx9SIia1R1fGO/Ybduhjonv5SwECE5Piqg5T4548mAlmeMadoPfvADfvCDHwBwxx13MHDgQAD69OlDdnY2/fr1Izs7m969ex8yXWlpKU8//TRvvvkmZ599NosXL+b555/nueee49prrz1k3KioKEJDQwHnoa1Ro0bx8ccfHxbL7NmzefXVVxkzZgwLFixg+fLlzcbumxhCQ0MpKyujtraWxMRE0tMbb8HYi9t1u/eloYJS+iTGEBoS2B92aNJQhiYNDWiZxpjG1V3y2b17N6+88gqXXXYZANOnT+fpp58GnCP/GTNmHDLdn/70J26++WbCw8MpKytDRAgJCam/rt+UESNGsH///vpEUFVVxYYNzp3vRUVF9OvXj6qqKp577rn6ac4880z+/ve/A84ZxMGDB5ssPyEhgbS0NP79738DTuJZu3YtAKeccgoLFy4EOKT8I9W9E4EHt44CvLP9Hd7Z/k7AyzXGHO7iiy9m5MiRXHjhhcybN4+kpCQAbrvtNt5++22GDx/O22+/zW23fdMSflZWFqtXr65PDrfccgsnn3wyTz/9NN/7XvPPvUZERLBo0SJ++ctfMmbMGMaOHVtf+XvvvfcyceJEpk6dyjHHHFM/zcMPP8z777/Pcccdx7hx4+oTR1Oee+45nnjiCcaMGcOoUaNYvHhxfTnz5s1jwoQJFBYWtv7HakK3rSNQVS7+81t8+7gB3HhuYJ8EtjoC0110hDoCczirI/BTUVkVJRXV9PfgjOCZ7zwT8DKNMcYr3TYRZLm3jvbvGRvwsgf1GNTySMYY00F02zqCrAMlQOCfIQB4I+MN3sh4I+DlGmOMF7rtGUF2/XsIAp8I7v/QeZR92rBpAS/bGGMCrdsmgj15xfTuEU1keGjAy144c2HAyzTGGK9020Swc18Rqb3jPSm7b1zgHv02xhivdcs6guqaWvbkFpPay5tEsHTL0vqGs4wxXdvy5cu54IILDuufnp7OsmXL2lTm73//+/ruQDRX3ZJumQi+PlBCda2S2iswbXk39ODHD/Lgxw96UrYxpvWqq6vbfZ7NJYKW4vFNBO2hW14a2pLlNER1VN8enpS/6JJFnpRrTIc3Zcrh/S65BK6/HkpL4bzzDh8+e7bzyc2FmTMPHdZCWz3gPM373HPPMWjQIFJSUhg3bhy33norU6ZMYfLkyXz00UdMnz6dmTNncs0117B//3569erFU089xeDBg5k9ezYXXHABM915x8XFUVxczPLly7n77rtJSUlh/fr1jBs3jmeffRYR4Y033mDu3LmkpKRw4oknHhZTZWUld911F2VlZXz44YfcfvvtbNq0iaysLHbu3ElKSgpnn302q1ev5tFHnXecX3DBBdx666288cYblJWVMXbsWEaNGsXvfve7+uaqV65cyYABA1i8eDHR0YFrKLNbnhGs23mAhOhwBnt0RpASk0JKTIonZRtjvrF69WpefvllvvjiC1555RUatjpQUFDAihUruOWWW7jxxhu56qqrWLduHZdffjk33XRTi+V/8cUXPPTQQ2zcuJHt27fz0UcfUV5ezrXXXsvSpUv54IMPyMk5/L0JERER3HPPPVx66aWkp6dz6aWXArBmzZr6xu2acv/99xMdHU16enp9e0JNNVcdKN3ujKBWlS925HLckGRCPHrpdt2LNi469iJPyjemw2ruCD4mpvnhKSl+nQH4+vDDD5kxY0b90fGFF154yPC6HTDAxx9/zCuvOP+bV155Jb/4xS9aLP+kk06qb8107Nix7Ny5k7i4ONLS0upfZHPFFVcwf/58v+KdPn16m47k09LSGDt2LADjxo1j586drS6jOZ6eEYjINBHZIiIZInJbI8MjReRFd/inIpLqZTwAX+zIJbeonG+N7OfZPB759BEe+fQRz8o3xjhaaistNrbplgPqmnMOCwujtra2vrzKysr6cRo2E113bb+tTUH7xuM7X4Dy8vImp2sqjkDxLBGISCgwDzgXGAlcJiIjG4z2AyBfVYcBfwH+6FU8ABVVNTz13haS4yOZdHQfz+azeNZiFs9a7Fn5xhjHqaeeytKlSykvL6e4uJj//ve/TY47efLkQ5pwPvXUUwFITU1lzZo1gPN6y6qqqmbnecwxx7Bjxw62bdsGwAsvvNDoePHx8RQVFTVZTmpqKunp6dTW1rJnzx4+++yz+mHh4eEtxhFIXp4RnARkqOp2Va0EFgIzGowzA3ja7V4EnClevHUBeP2L3Vz5yHtszS7k+nNGefIgWZ0eUT3oEeVNRbQx5hsTJkxg+vTpjBkzhosuuojx48fTo0fj/3uPPPIITz31FMcffzzPPPMMDz/8MADXXnstK1as4KSTTuLTTz9t9iwCnJfUzJ8/n/PPP59TTz2VIUOGNDreGWecwcaNGxk7diwvvvjiYcNPOeUU0tLSOO6447j11lsPqXSeM2cOxx9/PJdffrm/P8UR8awZahGZCUxT1R+6368EJqrqjT7jrHfHyXS/b3PHyW1Q1hxgDsDgwYPH7dq1q9XxfLZ1H8s3ZDF1zEBOSPO2Irfu5dyXjr60hTGN6dw6QjPUxcXFxMXFUVpayre+9S3mz5/f6J083UlHaoa6sSP7hlnHn3FQ1fnAfHDeR9CWYE4a3puThvduecQA+Ptq501ElgiM8d6cOXPYuHEj5eXlXH311d0+CbSFl4kgE/Btj3kgkNXEOJkiEgb0AA54GFO7WHZ5254mNMa0XnO3Yhr/eFlHsAoYLiJpIhIBzAKWNBhnCXC12z0TeE872yvTGhETHkNMeOBbNTWmI+oC/7JdSlvWh2eJQFWrgRuBN4FNwEuqukFE7hGR6e5oTwDJIpIB/Aw47BbTzujZdc/y7Lpngx2GMZ6LiooiLy/PkkEHoark5eURFRXVqum67TuLvWTvLDbdRVVVFZmZmc3eA2/aV1RUFAMHDiQ8PPyQ/vbO4nb29pVvBzsEY9pFeHg4aWlpwQ7DHCFLBB4IDw1veSRjjOkgumWjc15bkL6ABekLgh2GMcb4xRKBBywRGGM6k05XWSwi+4HWP1rsSAFyWxyra7Fl7h5smbuHI1nmIaraq7EBnS4RHAkRWd1UrXlXZcvcPdgydw9eLbNdGjLGmG7OEoExxnRz3S0R+Pcaoa7Flrl7sGXuHjxZ5m5VR2CMMeZw3e2MwBhjTAOWCIwxppvrNolARKaJyBYRyRCRTtvKqYgMEpH3RWSTiGwQkZvd/j1F5G0R2er+TXL7i4g84i73OhE50aesq93xt4rI1U3Ns6MQkVAR+UJEXnO/p4nIp278L7rNnSMike73DHd4qk8Zt7v9t4jIOcFZEv+ISKKILBKRze76ntTV17OI/NTdrteLyAsiEtXV1rOIPCki+9w3NNb1C9h6FZFxIvKlO80jIn68/ldVu/wHCAW2AUOBCGAtMDLYcbVxWfoBJ7rd8cBXwEjgT8Btbv/bgD+63ecBr+O8De5k4FO3f09gu/s3ye1OCvbytbDsPwOeB15zv78EzHK7/wFc53ZfD/zD7Z4FvOh2j3TXfSSQ5m4TocFermaW92ngh253BJDYldczMADYAUT7rN/ZXW09A98CTgTW+/QL2HoFPgMmudO8DpzbYkzB/lHa6YefBLzp8/124PZgxxWgZVsMTAW2AP3cfv2ALW73Y8BlPuNvcYdfBjzm0/+Q8TraB+cNd+8C3wZeczfyXCCs4TrGeQfGJLc7zB1PGq533/E62gdIcHeK0qB/l13PbiLY4+7cwtz1fE5XXM9AaoNEEJD16g7b7NP/kPGa+nSXS0N1G1idTLdfp+aeCp8AfAr0UdVsAPdv3Quam1r2zvabPAT8Aqh1vycDBeq8AAkOjb9+2dzhhe74nWmZhwL7gafcy2H/FJFYuvB6VtWvgQeA3UA2znpbQ9dez3UCtV4HuN0N+zeruySCxq6Rder7ZkUkDngZmKuqB5sbtZF+2kz/DkdELgD2qeoa396NjKotDOs0y4xzhHsi8HdVPQEoofk3+HX6ZXavi8/AuZzTH4gFzm1k1K60nlvS2mVs07J3l0SQCQzy+T4QyApSLEdMRMJxksBzqvqK23uviPRzh/cD9rn9m1r2zvSbnAJMF5GdwEKcy0MPAYkiUvdODd/465fNHd4DOEDnWuZMIFNVP3W/L8JJDF15PZ8F7FDV/apaBbwCTKZrr+c6gVqvmW53w/7N6i6JYBUw3L37IAKnYmlJkGNqE/cOgCeATar6fz6DlgB1dw5cjVN3UNf/Kvfug5OBQvfU803gbBFJco/Eznb7dTiqeruqDlTVVJx1956qXg68D8x0R2u4zHW/xUx3fHX7z3LvNkkDhuNUrHU4qpoD7BGREW6vM4GNdOH1jHNJ6GQRiXG387pl7rLr2UdA1qs7rEhETnZ/w6t8ympasCtN2rFy5jycO2y2Ab8KdjxHsByn4pzqrQPS3c95ONdG3wW2un97uuMLMM9d7i+B8T5lXQNkuJ/vB3vZ/Fz+KXxz19BQnH/wDODfQKTbP8r9nuEOH+oz/a/c32ILftxNEeRlHQusdtf1qzh3h3Tp9Qz8FtgMrAeewbnzp0utZ+AFnDqQKpwj+B8Ecr0C493fbxvwKA1uOGjsY01MGGNMN9ddLg0ZY4xpgiUCY4zp5iwRGGNMN2eJwBhjujlLBMYY081ZIjBdlogsFxHPX24uIje5rYM+16D/eBF5xO2eIiKTAzjPVBH5XmPzMqa1wloexZjuR0TC9Jv2bVpyPc696jt8e6rqapznAMB5/qEYWBmgGFKB7+G0xtpwXsa0ip0RmKByj2w3icjjbjv0b4lItDus/oheRFLcJiYQkdki8qqILBWRHSJyo4j8zG2c7RMR6ekziytEZKU47duf5E4f67YJv8qdZoZPuf8WkaXAW43E+jO3nPUiMtft9w+cB56WiMhPG4w/RURecxsHy7wepQAAAz9JREFU/DHwUxFJF5HTRKSXiLzsxrBKRE5xp7lbROaLyFvAv9zf5wMR+dz91J1V3A+c5pb307p5uWX0dH+fde7vcbxP2U+6v+t2EbnJ5/f4r4isdZft0iNbq6bTCfZTdvbp3h+cI9tqYKz7/SXgCrd7Oe6TlEAKsNPtno3zNGU80Aun1ckfu8P+gtMQX930j7vd38Jt9hf4vc88EnGeOI91y83EfaqzQZzjcJ7sjAXigA3ACe6wnUBKI9NM4ZunoO8GbvUZ9jxwqts9GKfJkLrx1vBNm/wxQJTbPRxY3bDsRub1V+A3bve3gXSfslfiPK2bAuQB4cDFdb+TO16PYG8X9mnfj10aMh3BDlVNd7vX4CSHlryvqkU47aoUAkvd/l8Cx/uM9wKAqv5PRBJEJBGnXZbpInKrO04Uzs4Y4G1VPdDI/E4F/qOqJQAi8gpwGvCFPwvYiLOAkfLNy6MSRCTe7V6iqmVudzjwqIiMBWqAo/0o+1ScnTuq+p6IJItID3fYf1W1AqgQkX1AH5zf7AER+SNOMvmgjctkOilLBKYjqPDprgGi3e5qvrl8GdXMNLU+32s5dLtu2IZKXVO9F6vqFt8BIjLx/7d39yoNREEYht8ROwt770HEn9LCSzCVhUhqQYXcgRjBJrVYWFh4B3YpLUQtRFJYWNmJIIgIaiGOxZzFENYkaxVyvqdLyBl208zuDMwQ457LDF73V80EsSzlo/vLlBi6r6EBPAFz6cznELH7jSLu/a8n3f3ezBaJmVUHZtZ2972h7kLGgnoEMsoeiJIM/E6frGoNwMyWicmNr8Tkxu00nREzmx8izjmwmiZjTgE1oMqT8xtRyiq0ga3iQ3riLzMNPLr7N7BBrF0ti9d7resp7grw7H12VpjZDPDu7qfEYpiFv34r40mJQEZZC9g0swuipv0fL+n8ETHlEaBJlFw6FgvEm4OCuPsNcEJMubwCjt29SlnoDKgVzWJgB1hKDd07oplc5hCom9klURYq3hY6wFdq8DZ6zuwWsYmmcp3+ZoFrM7slpnbuV7gvGQOaPioikjm9EYiIZE6JQEQkc0oEIiKZUyIQEcmcEoGISOaUCEREMqdEICKSuR/8koJQN/l0HQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "for i in range(rep):\n",
    "    plt.plot(mi_list[i,:],color='steelblue')\n",
    "    for t in range(mi_list[i].shape[0]):\n",
    "        if (mi_list[0,t]>.9*mi):\n",
    "            plt.axvline(t,label='90% reached',linestyle=':',color='green')\n",
    "            break\n",
    "plt.axhline(mi,label='ground truth',linestyle='--',color='red')\n",
    "plt.title(\"Plot of MI estimates against number of iteractions\")\n",
    "plt.xlabel(\"number of iterations\")\n",
    "plt.ylabel(\"MI estimate\")\n",
    "plt.legend()\n",
    "plt.savefig('./results/HGaussian_MINEE/MutualInfo.png')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
